目录
本文主要覆盖了肿瘤样本bulk RNA-seq数据常见的分析步骤,并从实践角度出发,较为具体地介绍了每一步骤依赖的工具和数据集。
另外,尽管本文适用于肿瘤样本,但其中的一些步骤,如基础分析、差异表达等是普适的,并不局限于肿瘤样本。
基础分析
1. 质控(reads)
- base call: 通常使用Q分数 (Phred quality score)衡量测序平台的准确性。
fastqc
, - adaptors:
Cutadapt
, - species/vector/ribosomal contamination:
- library complexity:
- mapping rates:
RSeQC
>bam_stat.py
- library stranded or not:
RSeQC
>infer_experiment.py
,
推荐工具:fastp
, RseqQC
Q分数 表示特定碱基被测序仪测错的概率。
2. 比对
推荐工具:STAR
首先,准备比对所需的基因组索引文件。索引文件获取方式以下两种:①使用STAR
生成,该操作依赖基因组序列和基因组注释两个数据集,推荐从Gencode网站下载。②从网络下载别人已经建好的索引文件(不推荐)。
如已有基因组索引文件,请跳过该步骤。
# 生成基因组索引
STAR \
--runMode genomeGenerate \ #运行index生成任务
--genomeDir /path/to/genomeDir \ #存储index路径
--genomeFastaFiles /path/to/genome/fasta1 ... \ #基因组参考序列文件
--sjdbGTFfile /path/to/annotations.gtf \ #转录本注释文件
--runThreadN 8 \ #指定使用线程数
--sjdbOverhang 150 #指定剪接点侧翼序列长度
然后,将样本测序数据比对到参考基因组。
STAR
--genomeDir /path/to/genomeDir \ #存储index路径
--readFilesIn /path/to/read1 ... \ #待比对序列文件,FASTQ或FASTA格式。支持多样本输入。
--readFilesCommand zcat \ #解压序列文件的明令
--outFileNamePrefix . #输出前缀
--sjdb* #version > 2.4.1a
--outFilterType BySJout #
--outFilterMultimapNmax 20 #
--outReadsUnmapped None
--alignSJoverhangMin 8 #
--alignSJDBoverhangMin 1 #
--outFilterMismatchNmax 999 #最大错配数
--outFilterMismatchNoverReadLmax 0.04 #最大错配比例
--outSAMstrandField intronMotif
--outSAMunmapped Within
--outSAMheaderHD @HD VN:1.4
--outSAMtype BAM Unsorted
--alignIntronMin 20 #内含子最小长度
--alignIntronMax 1000000 #内含子最大长度
--alignMatesGapMax 1000000 #mates间最大距离
--alignSJDBoverhangMin 10
--alignSJstitchMismatchNmax 5 -1 5 5
--runThreadN 8 \ #使用线程数
--chimSegmentMin 12 \
--chimJunctionOverhangMin 12 \
--chimNonchimScoreDropMin 10
--chimOutJunctionFormat 1 \
--peOverlapNbasesMin 12
--peOverlapMMp 0.1
--genomeLoad NoSharedMemory
--quantMode TranscriptomeSAM
3. 质控(alignment)
推荐工具:RseQC
是用于检查比对质量的常见工具,包括多个质量评估脚本,用户可以按需使用。
首先,为了减少运行时间,可以从BAM中选取部分比对进行质量评估。
# 1) 从比对采样50%
samtools view -s 0.5 -b sample.sorted.bam > sample_p50.bam
# 2) 选择管家基因的比对
bedtools intersect -a sample_p50.bam -b hk.bed > sample_p50_hk.bam
# 3) 建索引
samtools index sample_p50_hk.bam > sample_p50_hk.bam.bai
然后,使用RseQC
中相应的脚本评估质量。每个脚本依赖的数据集分别在对应命令下方列出。
bam_stat.py -i sample_p50_hk.bam #比对统计
geneBody_coverage.py -i sample_p50_hk.bam -r refGene.bed #reads在gene body覆盖分布
#1) 下载:https://sourceforge.net/projects/rseqc/files/BED/
read_distribution.py -i sample_p50_hk.bam -r refGene.bed #reads在基因组不同区域覆盖分布
RseQC质量指标:
- Transcript integrity number (TIN):整个基因组中具有均匀reads覆盖的转录本的百分比,用于衡量RNA完整性。
RseQC
会计算每个转录本的TIN,并报告样本中所有转录本TIN的平均数、中位数和标准差,中位数(medTIN)常用来表示每个样本的 RNA 完整性。 - Gene body coverage:显示了reads在基因体上的覆盖度。
- Read distribution:总结了比对到不同基因组区域(如外显子和内含子)的reads比例。
低质量样本的比对率低、完整性低或reads分布异常,应在下游分析中剔除。
4. 定量
定量主要包括计数和标准化两个步骤。标准化是从计数结果中去除测序深度、基因长度等非表达因素的影响,具体方法包括FPKM,TPM等。
下表总结了一些常见的用于定量的工具。
input | level | output | |
---|---|---|---|
featureCount | bam | 基因 | counts |
RSEM | 。 | 基因;转录本 | TPM, [RF]PKM, |
Salmon | 。 | 转录本 | TPM |
# 使用gtftools
HSA_GTF=Homo_sapiens.GRCh38.90.gtf
grep -v '#' $HSA_GTF | awk '$1=1' OFS='\t' > $HSA_GTF.tmp
gtftools.py -l Homo_sapiens.GRCh38.90.gtf.genelength Homo_sapiens.GRCh38.90.gtf.tmp
#
#gawk '$3~/exon/ && match($0, /gene_id "(.*)"; transcript_id "([^"]+)";/, a) {print $1,$4,$5,a[1],a[2]}' OFS='\t' $HSA_GTF |sort -k1,1V -k2,2n |bedtools merge -i - > tmp
Using GTF tools to get gene lengths
GTFtools
5. 样本合并
下游分析的输入一般需要将每个样本的分析结果合并,这通过简单脚本就可以方便的实现。
差异表达
1. 质控(cohort)
2. 差异分析
推荐工具:Deseq2
, EdgeR
, wilconxon rank-sum,
转录丰度文件,read计数矩阵,htseq-count文件,SummarizedExperiment对象都可以作为DESeq2
的输入文件。
# 示例1:Salmon输出 & tximport(用于将Salmon/Sailfish/kallisto/RSEM输出转换为转录丰度文件)
library("tximport")
library("readr")
library("tximportData")
library("DESeq2")
# 加载样本信息表
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
samples$condition <- factor(rep(c("A","B"),each=3)) #该列存储分组信息
#sample$condition <- relevel(sample$condition, ref="B") #设置B为比较的基线/野生型
rownames(samples) <- samples$run
#samples[,c("pop","center","run","condition")]
# 加载salmon结果
files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
# 将salmon结果转换成Deseq2能用的转录丰度文件
txi <- tximport(files, type="salmon", tx2gene=tx2gene)
# 用txi对象和样本信息构建DESeqDataSet对象
ddsTxi <- DESeqDataSetFromTximport(txi,
colData = samples,
design = ~ condition) #校正批次效应:design = ~ batch + condition,其中batch为样本信息表中记录批次信息的列名
dds <- DESeq(ddsTxi) #等效以下三步
#ddsTxi <- estimateSizeFactors(ddsTxi) # 测序深度归一化
#ddsTxi <- estimateDispersions(ddsTxi) # 基因离散度估计
#ddsTxi <- nbinomWaldTest(ddsTxi) # # 负二项GLM拟合 & Wald检验
DESeq2
输出值包括:
- baseMean:所有样本标准化后计数的均值;
- log2FoldChange:分组间倍数变化的log2转换值;
- lfcSE:log2FoldChange的标准误;
- stat:Wald统计量;
- pvalue:对Log2FoldChange估计的Wald检验的p值
- padj:FDR调整p值(BH方法)
3. 可视化(差异)
火山图的x轴代表log2FC,上调基因的log2FC为正,下调基因的log2FC为负。y轴代表 -log10(pvalue)。
富集分析
基因集富集分析(GSEA)可以检测预定义基因集(如DEG等)是否影响了特定通路、分子功能、细胞组分或过程。
推荐工具:clusterProfile
常见基因集:KEGG canonical pathways, GO(BP/CC/MF), MSigDB Hallmark gene sets, …
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(org.Hs.eg.db)
Gene Set Enrichment Analysis with ClusterProfiler
肿瘤免疫
1. 免疫组库
当肿瘤被T/B细胞浸润,如果其受体T/BCR能识别和结合肿瘤相关抗原,T/B细胞就会被激活。T/BCR具有高度多样性,以允许不同的T/B细胞识别不同的病原体或肿瘤相关抗原。T/BCR的多样性部分是由于V(D)J重组。
CDR3是T/BCR的一个高度可变区,是这些受体抗原结合部位的核心。免疫组库(描述肿瘤样本中 TCR/BCR的CDR3序列)对于量化T/B细胞多样性和克隆性至关重要。这些指标是重要的肿瘤免疫特征和免疫疗法反应的关键指标。
推荐工具:TRUST4
V(D)J重组会产生与参考基因组不一致的序列。TRUST
使用比对到参考基因组V/D/J区域的reads和无法比对到参考基因组的reads,推断从头组装的CDR3序列。
./run-trust4 \
-b example.bam \
-f hg38_bcrtcr.fa \ #包含V/D/J/C基因序列和坐标的fa文件
--ref human_IMGT+C.fa #来源于IMGT数据库的V/D/J/C基因参考fa文件
…
2. 免疫浸润
3. 免疫响应
4. 新抗原预测
肿瘤突变会导致蛋白质发生改变,这些蛋白质可能成为新抗原,并可能引起免疫反应。
适应性免疫依赖细胞表面呈现的peptide-MHC。人类 MHC 有两大类,即 I 类和 II 类(分别由 HLA-I 和 HLA-II 等位基因编码)参与抗原呈递。MHC I 类蛋白复合物在所有正常细胞上表达,呈递内部降解的蛋白质,而 MHC II 类蛋白复合物通常在APC上表达,呈递经过处理的外部抗原。MHC I 类蛋白复合物的抗原结合部分由 HLA-I 基因 HLA-A、HLA-B 和 HLA-C 编码。MHC II 类蛋白质的抗原结合部分由 HLA 复合物 DP、DQ 和 DR 区域中的 α 和 β 基因编码,包括 HLA-DQA1、HLA-DQB1 和 HLA-DRB1。目前基于比对的HLA分型方法需要输入 DNA 或 RNA 测序,只能预测 HLA-I 类或 HLA-II 两类。例如,arcasHLA
、HLAProfiler
和seq2HLA
就是通过RNA-seq数据进行高分辨率HLA分型。OptiType
和 PHLAT
是用于RNA、WES和WGS数据的HLA鉴定工具。HLA 参考序列可从 ImMunoGeneTics (IMGT) 数据库中获取。
新抗原预测算法通常使用NetMHCpan
、MHCflurry
和SMMalign
等工具整合肽/MHC 结合预测。
# 使用arcasHLA对RNA-seq数据进行HLA-I/II分型
## 提取reads
arcasHLA extract \
${sample}.sorted.bam \
-t 16 \
--sample $sample \
-o $out_np_dir
## 使用提取reads进行分型
arcasHLA genotype \
${sample}.extracted.1.fq.gz \
${sample}.extracted.2.fq.gz \
-g A,B,C,DQA1,DQB1,DRB1 \
-t 16 \
-o $out_np_dir
使用体细胞突变和 HLA 等位基因信息鉴定新抗原。体细胞突变检测的标准方法是分析肿瘤与配对正常组织的WES/WGS数据,常见的突变检测工具包括MuTect
, Varscan2
等。其中Varscan2
可以应用于具有良好覆盖率的 RNA-seq 数据,不过预计突变检测仍会比 DNA 测序产生的突变检测噪音大。之后,VEP
可用于注释基因、转录本和调控区的变异。
微生物组
在不同癌种中,肿瘤内微生物群的组成和丰度具有高度异质性。借助高通量测序技术,能够探索肿瘤宿主和宿主体内各种微生物的基因组,进而通过分析其与疾病表型、治疗效果等关系,阐明疾病发展的新机制,为发现治疗目标提供新见解。
一般使用16S rRNA和鸟枪法测序完成微生物系统发育和分类。研究人员还开发了使用基于参考基因组的高通量测序方法,如PathSeq
和Centrifuge
。PathSeq
利用 BWA-MEM 将非宿主reads与预先定义的微生物序列进行比对,而Centrifuge
则先建立基因组序列的压缩索引。虽然两者都能实现对微生物群的分类和丰度估计,但在内存使用和速度方面,Centrifuge
优于PathSeq
。
Centrifuge
使用FM-index来存储数据库中的所有序列信息。Centrifuge
会在reads和参考数据库之间贪婪寻找长匹配序列,并根据匹配长度进行打分。得分最高的种或属将是reads的分类结果。reads可以被分给多个物种,Centrifuge
使用EM算法来量化已识别物种的丰度。由于微生物组参考数据库的规模不断扩大,Centrifuge
使用同属物种之间的高序列相似性,可以压缩参考数据库。
# 使用RNA-seq进行微生物分类
centrifuge \
-x ref_files/centrifuge_index/p_compressed+h+v \ #参考索引下载:https://ccb.jhu.edu/software/centrifuge/manual.shtml#nt-database
-p 16 \
--host-taxids 9606 \
-1 ${sample}_1.fastq.gz \
-2 ${sample}_2.fastq.gz \
-s ${sample}_classification.txt.gz \
--report-file ${sample}_report.txt
结果输出:
- name:微生物名
- taxID:物种ID
- taxRank:分类单位
- genomeSize:
- numReads:
- numUniqueReads:
- abundance:
参考
Yang Liu, Jennifer Altreuter, Yang Lin. Tutorial of RNA-seq tumor immunity analysis
RNA-seq免疫分析流程RIMA的教程