使用工具
fastp(質(zhì)控), hisat2(比對), samtools(sam文件轉(zhuǎn)bam文件), featureCounts(count計數(shù)), DESeq2(差異分析)
環(huán)境配置
使用conda配置環(huán)境, 安裝fastp, hisat2, samtools, subread。featureCounts整合在subread中。DESeq2為R包。使用DESeq2進(jìn)行差異分析前的流程都在linux環(huán)境下進(jìn)行，DESeq2在R環(huán)境下運行。
分析流程
- 拿到數(shù)據(jù)后首先需要對數(shù)據(jù)進(jìn)行質(zhì)控，這里使用fastp。這篇文章介紹的比較清楚：[鏈接1]。
- 使用hisat2進(jìn)行序列比對，輸入文件為測序數(shù)據(jù)fastq(gz)文件，輸出文件為sam文件，額外需要基因組索引文件，需要到hisat2官網(wǎng)[鏈接2]下載。
hisat2 -p 8 -x ref/grch38/genome -1 sample_1.fq.gz -2 sample_2.fq.gz -S sample.sam --summary-file sample.hisat2.summary
# -p: 線程數(shù)
# -x: 基因組索引前綴。所下的基因組索引為多個文件，索引前綴到genome為止。
# -1/-2: fastq輸入文件。當(dāng)輸入為單端測序時使用-U 指定輸入。單端或雙端的輸入都可使用逗號分隔輸入多個樣本，或使用多次-1 -2 / -U 指定多個輸入。e.g.: -U sample1.fq.gz,sample2.fq.gz -U sample3.fq.gz
# -S: 輸出sam文件路徑。
# --summary-file: 生成summary文件。
- 轉(zhuǎn)換成bam文件
sam文件為通用的比對數(shù)據(jù)存儲格式，記錄了reads的比對信息，其內(nèi)容為純文本，因此大小通常十分大。為解決這個問題，sam文件的壓縮格式bam文件被設(shè)計了出來，使文件大小被大大壓縮。這里通過samtools進(jìn)行bam文件轉(zhuǎn)換：
samtools sort -@ 8 -o sample.bam smaple.sam
# sort: 進(jìn)行排序
# -@: 線程數(shù)
# -o: 輸出bam文件
# 最后一項為輸入文件
- 使用featureCounts進(jìn)行計數(shù)，需要基因組注釋文件，可在gencode網(wǎng)站下載[鏈接3]。
featureCounts -p -T 8 -a ref/annotation.gtf.gz -o featurecounts.out sample1.bam sample2.bam ...
# -p: 此參數(shù)雙端測序時使用
# -T: 線程數(shù)
# -a: 基因組注釋文件
# -o: 輸出文件
# 最后為bam文件，可指定多個輸入
- 使用DESeq2進(jìn)行差異分析
差異分析在R環(huán)境中進(jìn)行。
1) 構(gòu)建dds對象
dds <- DESeqDataSetFromMatrix(countData, colData, design, tidy = False, ignoreRank = False)
# countData: 非負(fù)整數(shù)count矩陣，count數(shù)據(jù)從第一列開始，列名為樣本名，行名可指定為gene id。
# colData: 行名為countData的列名且至少1列的DataFrame，每一列為樣本分組信息和其他批次信息。
# design: 指定colData中的分組變量，e.g.: design = ~batch+type (batch和type均為colData的列名)。DESeq2將默認(rèn)使用最后一個分組變量進(jìn)行差異分析。
# tidy: 邏輯變量，countData第一列是否為行名。
# ignoreRank: use of this argument is reserved for DEXSeq developers only. Users will immediately encounter an error upon trying to estimate dispersion using a design with a model matrix which is not full rank.
Note: countData列名和colData行名必須相同。
2) 差異分析---differential experssion analysis
dds_deseq <- DESeq(dds)
這一步將經(jīng)過3步分析:
1. estimation of size factors
2. estimation of dispersion
3. Negative Binomial GLM fitting and Wald statistics
3) 提取結(jié)果
res <- results(dds_deseq, contrast = c("type","treatment","untreatment"))
# contrast: 3個元素的向量，第一個元素: colData中的分組列名，第二個元素: 計算log2FC時的分子項，第三個元素: 計算log2FC時的分母項。
res <- res[order(res$padj),]
# 根據(jù)padj排序
summary(res)
# 查看summary
?
4) DEseq標(biāo)準(zhǔn)化數(shù)據(jù)
除了進(jìn)行差異分析外，有時我們還需要直接得到標(biāo)準(zhǔn)化的表達(dá)矩陣，現(xiàn)有的標(biāo)準(zhǔn)化方法有CPM、TPM、RPKM/FPKM、TMM、DESeq標(biāo)準(zhǔn)化等。這篇文章[鏈接4]詳細(xì)介紹了各標(biāo)準(zhǔn)化的適用范圍以及DESeq標(biāo)準(zhǔn)化的計算方法。
代碼：
#創(chuàng)建dds對象
dds <- DESeqDataSetFromMatrix(countData, colData, design)
#計算量化因子
dds_SF <- estimateSizeFactors(dds)
#進(jìn)行標(biāo)準(zhǔn)化
normalized_counts <- counts(dds_SF, normalized=TRUE)
另外DESeq提供了其他兩種標(biāo)準(zhǔn)方法VST和rlog，但目前對這兩種標(biāo)準(zhǔn)化方法的理解還不夠深入。
兩種方法的差異及如何選擇：
For genes with high counts, both the VST and the rlog will give similar result to the ordinary log2 transformation of normalized counts. For genes with lower counts, however, the values are shrunken towards a middle value. The VST or rlog-transformed data then become approximately homoskedastic (more flat trend in the meanSdPlot), and can be used directly for computing distances between samples, making PCA plots, or as input to downstream methods which perform best with homoskedastic data.
Which transformation to choose?
The VST is much faster to compute and is less sensitive to high count outliers than the rlog. The rlog tends to work well on small datasets (n < 30), potentially outperforming the VST when there is a wide range of sequencing depth across samples (an order of magnitude difference). We therefore recommend the VST for medium-to-large datasets (n > 30).
具體介紹見這里[鏈接5]
兩種轉(zhuǎn)換方法可分別用vst和rlog函數(shù)實現(xiàn)：
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
