Bulk RNA-Seq 差异表达分析流程
以前写过不少零散的 RNA-Seq 分析文章，现在整理为流程，同时修改一些错误。
流程包含质控、比对、定量、差异分析。
流程概况
前处理
拿到原始 fastq 数据先进行前处理。前处理包含质控、比对和定量。质控采用 fastqc/fastp; 比对用 hisat2 或者不比对，用 salmon 直接定量；比对后用 featureCounts 进行 reads 定量，用 TPMCalculator 进行 TPM 定量。
fastp 质控
推荐 fastp 质控原始 fastq 数据，因为集成程度高包含了各项过滤同时速度很快。
fastp -i ${RawDir}/${Sample}_R1.fq.gz -o ${CleanDir}/${Sample}_R1.fq.gz \
-I ${RawDir}/${Sample}_R2.fq.gz -O ${CleanDir}/${Sample}_R2.fq.gz \
--compression 6 --report_title ${Sample} --json ${CleanDir}/${Sample}_fastp.json \
--html ${CleanDir}/${Sample}_fastp.html --detect_adapter_for_pe
一般来说 RNA-Seq 重复率会比较高。ATCG 碱基占比容易出现前面碱基不稳定。多查看几个样本表现，如果波动模式差不多说明是随机引物不随机导致的，可以选择移除或不移除。
ATCG 占比前面碱基有波动
hisat2 比对
比对推荐 hisat2 不推荐 bowtie2. 因为 bowtie2 对剪切不友好，即使提高 --maxins 参数允许更长的片段，比对率依旧比较低。
hisat2 --new-summary -k 4 --threads 4 --summary-file ${BamDir}/${Sample}_HISAT2.txt \
-x ${GRCh38} -1 ${CleanDir}/${Sample}_R1.fq.gz -2 ${CleanDir}/${Sample}_R2.fq.gz \
-S ${BamDir}/${Sample}.sam
# samtools 过滤
samtools view --threads 4 -F 256 -b ${BamDir}/${Sample}.sam | samtools sort \
--threads 4 -o ${BamDir}/${Sample}.bam -O BAM
samtools index ${BamDir}/${Sample}.bam
ENCODE 2016 指南认为 Long/Small RNA-Seq 文库应有 30M 以上比对 Reads/Mates.
hisat2 默认 -k 参数为 5(linear index) 或 10(graph index)，此时 MAPQ 数值将为 0 或 1 或 60. 60 为唯一比对（primary），0 为不比对，1 为多比对。
$ samtools view Test.bam | awk '{print $5}' | sort | uniq
0
1
60
在 samtools view 命令用 -F 256 移除次要比对，保留唯一比对。
$ samtools flag 256
0x100 256 SECONDARY
featureCounts 定量
featureCounts 定量基础单元叫 feature (如外显子 exon)，以 GTF 注释文件为例，其第三列为 feature 类型。多个 feature 可联合形成 meta-feature (如基因 gene) 并计算总聚的 reads 数目。由 -g 参数设置 meta-feature 对应哪个属性，比如 GTF 格式默认用 "gene_id" 标记。
featureCounts 默认不计算重叠于多个 (meta-)feature 的 reads，除非重叠的 feature 属于同一个 meta-feature 那么只计算一次。重叠于多个 meta-feature 对 RNA-Seq 实验建议不计算，因为无法得知该 reads 究竟来自哪个基因；对于 CHIP-Seq 建议计算。因此用 featureCounts 定量转录本是不合适的，因为同一基因的转录本有大量重合区域。
featureCounts 用 BAM 文件 "NH" 标签获取多比对信息。默认只计算唯一比对（primary），不计算多比对。用 -M 参数控制计算多比对；用 --primary 参数控制只计算唯一比对。
featureCounts 默认 reads 有一个碱基与 feature 重合就计算一次，可通过 --minOverlap, --fracOverlap 等参数设定阈值。同时 reads 上的任意 gap(insertions, deletions, exon-exon junctions or structural variants) 也算数。
featureCounts 输入的 BAM 文件不要求排序。下面 -p 和 --countReadPairs 参数指定是双端测序且统计 fragment 而不是 reads.
featureCounts -p --countReadPairs --primary -F GTF -t exon -g gene_id \
-a ${Annot} -o ${CountDir}/${Sample}.tsv ${BamDir}/${Sample}.bam
salmon 定量
salmon 是转录本定量软件，使用转录本定量主要优点一是更准确，比如样本同一基因使用不同转录本此时基因长度并不相等，直接基因定量不准确；二是方便进行可变剪切分析。
salmon 可以从 fastq 直接定量也可以从比对好的 bam 定量，本文使用 fastq 定量。
第一次使用 salmon 前先建立索引，建索引前先准备 decoys 文件。decoys 指不属于人基因组但是测序又往往会检测到的序列，有些 decoys 序列跟一些转录本很相似，明确 decoys 序列有助于准确定量。
用 GitHub - COMBINE-lab/SalmonTools: Useful tools for working with Salmon output 仓库的脚本可以制作 decoys 文件。
bash generateDecoyTranscriptome.sh -j 8 -a gencode.v35.annotation.gtf \
-g GRCh38.primary_assembly.genome.fa -t gencode.v35.transcripts.fa -o hsa_decoy
用 decoys 目录里的 fa 文件建立索引。
salmon index -t hsa_decoy/gentrome.fa -i hsa_transcripts_index \
--decoys hsa_decoy/decoys.txt -k 31 --gencode
设置 -k 31 适合 75 bp 及更长 reads，如果 reads 更短需要设置小点。或者比对率偏低时，可以将 -k 设小些，提高灵敏度。
定量前注意两点，一是 salmon 要求 reads 无序；二是注意文库类型。
ReadLibraryIllustration
如上图所示，根据 reads 方向将 RNA-Seq 文库分为不同类型。salmon LIBTYPE 参数包含三个部分，一是 reads 的相对方向，分别用 I/O/M 表示；二是文库有无链特异性，分别用 S/U 表示；三是，假如文库有链特异性，指明链方向，分别用 F/R 表示。如果第二部分是 U 就不需要设置第三部分。下面是这些字母代表的全称。
I = inward
O = outward
M = matching
S = stranded
U = unstranded
F = read 1 (or single-end read) comes from the forward strand
R = read 1 (or single-end read) comes from the reverse strand
A = automatically determine
salmon 用 F/R 表示链方向，其他软件可能用 F2R1/F1R2 形式。
下图是与 TopHat 比较
LYT_Salmon_TopHat
注意命令 -l 参数放 -1/-2 之前。
salmon quant -l IU -i ${IndexPath} -1 ${CleanDir}/${Sample}_R1.fq.gz \
-2 ${CleanDir}/${Sample}_R2.fq.gz -p 4 -g ${GRCh38GTF} \
-o ${SalmonDir}/${Sample}
转录本定量保存在 quant.sf 文件。有 -g 参数时会输出基因水平定量。
Name Length EffectiveLength TPM NumReads
ENST00000456328.2 1657 1372.105 0.034568 1.000
ENST00000450305.2 632 348.480 0.000000 0.000
ENST00000488147.1 1351 1066.105 3.604749 81.025
注意 NumReads 不一定是整数，是考虑唯一比对和多比对后对转录本相对丰度的估计。
TPMCalculator 定量
TPMCalculator 计算基因、转录本、外显子和内含子的 TPM 定量。其 TPM 定量公式如下， 是比对到
feature 片段数目；
是 feature 长度。
基因的转录本及外显子内含子有着复杂的关系，如下图所示软件进行了两种转换。第一种转换，根据外显子重合，创建包含所有外显子的基因模型；第二种转换创建“纯粹”的内含子区域，不被任何外显子覆盖。这些不重合的内含子和 feature 用来计算 TPM.
Gene_model
参数 -c 推荐设置为 reads 长度；参数 -p 设置只统计“正确”比对的双端测序数据。
因为比对 hisat2 设置了 -k 参数，因此 MAPQ 值只有 0/1/60 三种，用 -q 参数过滤比对质量 60 因此只保留 primary alignment. 如果前面经过 samtools 过滤了，这里不再需要过滤。
TPMCalculator -g ${GRCh38GTF} -b ${BamDir}/${Sample}.bam \
-c 150 -k gene_id -t transcript_id -q 60 -p -a
软件将输出结果至当前目录，包含 3 文件。
- *_genes.out | Include calculated values per Gene
- *_genes.uni | Include calculated values per non-overlapped features per Gene
- *_genes.ent | Include calculated values for all features per Gene
仔细检查基因输出结果可能会发现有的基因重复了，这是因为软件认为有重复的基因。但提供的 GTF 这个基因并没有重复，但是它存在不重叠的转录本，软件就认定为重复的基因，这算是软件的一个 bug 了，希望哪天能修复了。
$ grep "ENSG00000235538" kLEC96h-2_genes.out | awk '{NF=7;print}'
ENSG00000235538.3#1 chr6 163671576 163798848 127273 4 0.00366198
ENSG00000235538.3#2 chr6 163927121 164009979 82859 2 0.00281244
ENSG00000235538.3#3 chr6 164228330 164231609 3280 0 0
差异表达分析
差异分析用 DESeq2 包，定量采用 salmon 结果，需要 tximport 将 salmon 导入给 DESeq2. 如果用 featureCounts 定量，将 counts 整理成矩阵直接导入 DESeq2 分析。
tximport 将 salmon 定量的转录本表达计算成基因表达。这需要基因和转录本映射表（TxDb），TxDb 用 GenomicFeatures 包从 GTF 注释文件创建，然后保存到本地，以后就读取使用。
library(AnnotationDbi)
library(GenomicFeatures)
gtfPath <- file.path("/example_dir/gencode.v35.annotation.gtf")
metaName <- c("Source", "Version", "Species")
metaValue <- c("GENCODE", "v35", "Homo sapiens")
extraInfo <- data.frame(name=metaName, value=metaValue)
txd <- makeTxDbFromGFF(gtfPath, format = "gtf", dataSource = "gencode.v35.annotation.gtf", organism = "Homo sapiens", metadata = extraInfo)
# 保存 TxDb 对象到本地
txdbPath <- file.path("/example_dir/gencode.v35.annotation.TxDb.sqlite")
saveDb(txd, file = txdbPath)
用 makeTxDbFromGFF 读取 GTF 注释到 TxDb. 参数 dataSource 记载数据来源，可填可不填；参数 metadata 提供元信息，格式是两列的数据框，要求列名分别为 "name", "value".
> keytypes(txd)
[1] "CDSID" "CDSNAME" "EXONID" "EXONNAME" "GENEID" "TXID" "TXNAME"
> txToGene <- AnnotationDbi::select(txd, keys=keys(txd, "TXNAME"), keytype="TXNAME", columns=c("TXNAME", "GENEID"))
'select()' returned 1:1 mapping between keys and columns
> head(txToGene)
TXNAME GENEID
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000473358.1 ENSG00000243485.5
4 ENST00000469289.1 ENSG00000243485.5
5 ENST00000607096.1 ENSG00000284332.1
6 ENST00000606857.1 ENSG00000268020.3
txToGenePath <- file.path("/example_dir/gencode.v35.annotation.TxToGene.csv")
write.csv(txToGene, file = txToGenePath, quote = FALSE, row.names = FALSE)
选择转录本和基因两列，得到他们映射关系，并保存到本地 csv 文件，下次直接读取使用。注意第一列为转录本 ID 第二列为基因 ID 固定顺序，列名不作要求。
有 TxDb 和 salmon 定量文件和样品分组信息，就可以将数据导入到 DESeq2 分析。
导入必要的 R 包和文件路径。biomaRt 用于基因名转换，因为差异分析我习惯用 Ensembl ID 下游分析又常用 HGNC SYMBOL 或 ENTREZ ID.
library(tidyverse)
libraty(biomaRt)
library(tximport)
library(DESeq2)
groupPath <- file.path("/example_dir/SampleGroup.csv")
txGenePath <- file.path("/example_dir/gencode.v35.annotation.TxToGene.csv")
salmonDir <- file.path("/example_dir/Salmon")
ensembl <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
tximport 读取 salmon 定量，并保存一份 TPM 表达矩阵，以备下游分析需要。要注意给 salmon 文件路径命名，保证读取后矩阵的列跟文件是对应的。
sampleGroup <- read.csv(groupPath, header = TRUE, quote = "", row.names = 1)
sampleList <- rownames(sampleGroup)
sampleFiles <- file.path(salmonDir, sampleList, "quant.sf")
# 注意给路径命名
# 名字是表达矩阵样品名
names(sampleFiles) <- sampleList
txToGene <- read.csv(txGenePath, header = TRUE, quote = "", stringsAsFactors = FALSE)
txi <- tximport(sampleFiles, type = "salmon", tx2gene = txToGene)
# 理论上 TPM 不应该进行过滤一些基因，除非在所有样本 TPM 都为 0
tpm1 <- txi$abundance
tpm2 <- tpm1[rowSums(tpm1) > 0, ]
# 去除 Ensembl ID 版本号
short_id <- function(long_id) {
id_version <- strsplit(long_id, split = ".", fixed = TRUE) %>% unlist()
ensembl_id <- id_version[1]
return(ensembl_id)
}
# 取得所有表达基因，创建基因 ID 映射表
ensemblGenes <- sapply(rownames(tpm2), short_id, USE.NAMES = FALSE)
geneMap <- getBM(mart = ensembl, attributes = c("ensembl_gene_id", "entrezgene_id", "hgnc_symbol"),
values = ensemblGenes, filters = "ensembl_gene_id") %>%
as_tibble() %>%
arrange(entrezgene_id, desc(hgnc_symbol)) %>%
distinct(ensembl_gene_id, .keep_all = TRUE)
# 保存 TPM 矩阵到文件
tpm3 <- as_tibble(tpm2, rownames="gene_id") %>%
tidyr::separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
dplyr::select(-ensembl_gene_version) %>% dplyr::left_join(geneMap, by="ensembl_gene_id") %>%
dplyr::select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
tpmPath <- file.path("/example_dir/TPM.csv")
write_csv(tpm3, tpmPath)
导入到 DESeq2 并过滤低表达基因。过滤没有统一标准，我习惯要求至少在 n 样本 counts 不小于 x. 如果数据有不同批次，一般在 design 公式将批次列（因子）放前面，注意 DESeq2 不会进行批次效应移除，但分析时会区分哪些差异由批次效应引起，哪些由实验条件引起。
dds1 <- DESeqDataSetFromTximport(txi, colData = sampleGroup, design = ~ Group)
# 过滤低表达基因
# 一共 6 样品，要求至少 2 样品 read counts 大于等于 5
keepGene <- rowSums(counts(dds1) >= 5) >= 2
dds2 <- dds1[keepGene, ]
dds3 <- DESeq(dds2)
用 cook's 距离检测组内样本一致性。一般 RNA-Seq 测序至少有三个生物学重复，如果重复少于三个无法计算 cook's 距离。DESeq2 不支持技术重复，如果做了技术重复，用 collapseReplicates 函数合并。
ddsAssay <- assays(dds3)
cooks <- ddsAssay$cooks
boxplot(log10(cooks))
取得差异分析结果前保存表达矩阵，以备数据质控和下游分析。用 counts 函数取得原始的 counts 或标准化的 counts. 后面差异分析结果的 baseMean 列就是标准化的 counts 计算得到。函数 rlog 或 vst 将 counts 根据文库大小或其他标准化因子进行了标准化和转换到对数（log2）尺度，并控制方差独立于表达均值。参数 blind 控制转换过程是否考虑实验设计——即前面的 design 参数。设置 blind = FALSE 会考虑 counts 差异是否是实验条件导致，从而保留应有的差异避免过度压缩，得到的矩阵适合进行下游分析。
readCounts <- DESeq2::counts(dds3, normalized = FALSE) %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
normalizedCounts <- DESeq2::counts(dds3, normalized = TRUE) %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
# 适合质控
rLog1 <- rlog(dds3, blind = TRUE) %>% assay() %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
# 适合下游分析
rLog2 <- rlog(dds3, blind = FALSE) %>% assay() %>%
as_tibble(rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by = "ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
取得差异分析结果并保存 MA 图。其中 results 函数取得“原始”的差异表达数据，lfcShrink 将差异倍数“压缩”，两个数据 p 值相同——即统计检验的结论是相同的。“压缩”后的差异表达数据适合用于排序或可视化。要注意如果实验设计 design 包含了交互作用项，不要进行 lfcShrink 因为容易过度“压缩”。
参数 contrast 接受多种指定方式，下面代码是常见的两种方式。第一种是三个元素的向量，向量第一个元素是因子名字，第二个元素是差异分析分子项的水平，即常指的实验组（比较组），第三个元素是差异分析分母项的水平。举例来说，c("Condition", "B", "A") 表示条件 B 样本对比条件 A 样本的差异结果。第二种从 resultsNames() 结果选取一个字符串，如 "Group_48H_vs_0H" 表示 Group 因子项的 48H 样本比较 0H 样本结果。
rawRes <- results(dds3, contrast = c("Group", "48H", "0H"))
lfcRes <- lfcShrink(dds3, coef = "Group_48H_vs_0H", type = "apeglm", res = rawRes)
summary(rawRes)
# MA Plot
plotMA(rawRes, ylim = c(-5, 5), main = "Raw")
plotMA(lfcRes, ylim = c(-5, 5), main = "Shrink")
比较 MA 图和火山图，可以看出 lfcShrink 对低表达基因“压缩”影响更大。
MA_plot
Volcano
将差异分析结果转换为 tibble 并保存，下游分析根据需要进行过滤。一般来说取 "abs(log2FoldChange) >= 1, padj < 0.05" 作为差异表达基因。
rawDEGs <- as_tibble(rawRes, rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by="ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
lfcDEGs <- as_tibble(lfcRes, rownames = "gene_id") %>%
separate(col = gene_id, into = c("ensembl_gene_id", "ensembl_gene_version"), sep="\\.", remove = TRUE, extra = "drop", fill = "right") %>%
select(-ensembl_gene_version) %>%
left_join(geneMap, by="ensembl_gene_id") %>%
select(ensembl_gene_id, entrezgene_id, hgnc_symbol, everything())
差异基因结果基因 p 值可能是 NA. 包含下列情况：
- 该基因所有样本 counts 数都为 0, 那么 baseMean 列为 0, 此时 log2FC, p, padj 为 NA
- 该基因包含离群的 counts 数，此时 p, padj 为 NA
- 该基因 baseMean 非常低，此时 padj 为 NA.
相互作用项
上面提到实验设计包含相互作用项时不建议 lfcShrink 处理，相互作用项在官网 Interactions 有详细说明。
简单来说，如果认为因子间有相互作用，不是独立的。应在 design 参数添加相互作用项。
假设实验设计有 2 因子 genotype 值分别为 I/II/III 和 condition 值分别为 A/B.
如果 design 为 ~ genotype + condition 那么结果将是总体的 condition 导致的差异基因，排除了 genotype 因子的效应；如果是 ~ genotype + condition + genotype:condition 或者写为 ~ genotype * condition ，那么 condition 因子效应将拆分为主要的效应——即 condition 因子在 genotype 为对照组（I）的情况下的效应，及相互作用项genotypeII.conditionB 和 genotypeIII.conditionB 表示对比于 genotypeI 在 genotypeII/genotypeIII 条件下的 condition 效应的差异。
下面举个例子。基因 1 在三种 genotype 下不同 condition 导致的差异倍数几乎相同，那么相互作用项 genotype:condition 将接近于 0. 基因 2 在三种 genotype 下有不同的 condition 效应，那么 genotype 为 II/III 的 condition 效应将是主要的 condition 效应加上对应 genotype 的相互作用项。
Gene12
参考资料
Vera Alvarez, Roberto, et al. "TPMCalculator: one-step software to quantify mRNA abundance of genomic features." Bioinformatics 35.11 (2019): 1960-1962.
QC Fail Sequencing » MAPQ values are really useful but their implementation is a mess
Analyzing RNA-seq data with DESeq2
Rsubread/Subread Users Guide
Overview – Salmon: Fast, accurate and bias-aware transcript quantification from RNA-seq data
Bowtie 2: fast and sensitive read alignment
The decoy genome
https://www.biostars.org/p/456231/
Decoy Sequences V/S Target Sequences
Importing transcript abundance with tximport
