使用工具
fastp(质控), hisat2(比对), samtools(sam文件转bam文件), featureCounts(count计数), DESeq2(差异分析)
环境配置
使用conda配置环境, 安装fastp, hisat2, samtools, subread。featureCounts整合在subread中。DESeq2为R包。使用DESeq2进行差异分析前的流程都在linux环境下进行，DESeq2在R环境下运行。
分析流程
- 拿到数据后首先需要对数据进行质控，这里使用fastp。这篇文章介绍的比较清楚：[链接1]。
- 使用hisat2进行序列比对，输入文件为测序数据fastq(gz)文件，输出文件为sam文件，额外需要基因组索引文件，需要到hisat2官网[链接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: 线程数
# -x: 基因组索引前缀。所下的基因组索引为多个文件，索引前缀到genome为止。
# -1/-2: fastq输入文件。当输入为单端测序时使用-U 指定输入。单端或双端的输入都可使用逗号分隔输入多个样本，或使用多次-1 -2 / -U 指定多个输入。e.g.: -U sample1.fq.gz,sample2.fq.gz -U sample3.fq.gz
# -S: 输出sam文件路径。
# --summary-file: 生成summary文件。
- 转换成bam文件
sam文件为通用的比对数据存储格式，记录了reads的比对信息，其内容为纯文本，因此大小通常十分大。为解决这个问题，sam文件的压缩格式bam文件被设计了出来，使文件大小被大大压缩。这里通过samtools进行bam文件转换：
samtools sort -@ 8 -o sample.bam smaple.sam
# sort: 进行排序
# -@: 线程数
# -o: 输出bam文件
# 最后一项为输入文件
- 使用featureCounts进行计数，需要基因组注释文件，可在gencode网站下载[链接3]。
featureCounts -p -T 8 -a ref/annotation.gtf.gz -o featurecounts.out sample1.bam sample2.bam ...
# -p: 此参数双端测序时使用
# -T: 线程数
# -a: 基因组注释文件
# -o: 输出文件
# 最后为bam文件，可指定多个输入
- 使用DESeq2进行差异分析
差异分析在R环境中进行。
1) 构建dds对象
dds <- DESeqDataSetFromMatrix(countData, colData, design, tidy = False, ignoreRank = False)
# countData: 非负整数count矩阵，count数据从第一列开始，列名为样本名，行名可指定为gene id。
# colData: 行名为countData的列名且至少1列的DataFrame，每一列为样本分组信息和其他批次信息。
# design: 指定colData中的分组变量，e.g.: design = ~batch+type (batch和type均为colData的列名)。DESeq2将默认使用最后一个分组变量进行差异分析。
# 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)
这一步将经过3步分析:
1. estimation of size factors
2. estimation of dispersion
3. Negative Binomial GLM fitting and Wald statistics
3) 提取结果
res <- results(dds_deseq, contrast = c("type","treatment","untreatment"))
# contrast: 3个元素的向量，第一个元素: colData中的分组列名，第二个元素: 计算log2FC时的分子项，第三个元素: 计算log2FC时的分母项。
res <- res[order(res$padj),]
# 根据padj排序
summary(res)
# 查看summary
4) DEseq标准化数据
除了进行差异分析外，有时我们还需要直接得到标准化的表达矩阵，现有的标准化方法有CPM、TPM、RPKM/FPKM、TMM、DESeq标准化等。这篇文章[链接4]详细介绍了各标准化的适用范围以及DESeq标准化的计算方法。
代码：
#创建dds对象
dds <- DESeqDataSetFromMatrix(countData, colData, design)
#计算量化因子
dds_SF <- estimateSizeFactors(dds)
#进行标准化
normalized_counts <- counts(dds_SF, normalized=TRUE)
另外DESeq提供了其他两种标准方法VST和rlog，但目前对这两种标准化方法的理解还不够深入。
两种方法的差异及如何选择：
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]
两种转换方法可分别用vst和rlog函数实现：
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
