今天分享的学习笔记是一套转录组分析简单流程，适用于初学者入门阅读，从原始测序数据开始，经过质控、序列比对、定量表达、差异表达、功能富集等一系列分析步骤，最终获得基因表达信息，制作出火山图和功能富集图。本文所有数据都经过特殊修改，仅供学习参考使用。
”
转录组是在特定时空条件下细胞中基因转录表达产物，广义的转录组包括信使RNA，核糖体RNA，转运RNA及非编码RNA，狭义上是指所有mRNA的集合，转录组分析能够获得不同基因的表达情况。
1. 数据来源
假设有两个不同组织（PR和SR），每个组织各区三个样本，一共六个样本，利用illumina平台进行转录组测序，得到双端测序数据。数据原始格式为.fq,共有12条测序数据文件（每个样本产生两条）
PR1_2.fq PR2_2.fq PR3_2.fq SR1_2.fq SR2_2.fq SR3_2.fq PR2_1.fq PR3_1.fq SR1_1.fq SR2_1.fq SR3_1.fq
创建工作文件夹，并新建子步骤目录：00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis
mkdir RNA-Seq_Practice cd RNA-Seq_Practice mkdir 00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis #批量创建工作文件夹
2. 测序数据质量评估
利用fastQC软件对获得的fastq序列文件进行质量分析，生成html格式的结果报告，其中含有各项指标，以下用PR1样本为例。
fastqc *.fq #批量对fq后缀文件运行fastqc程序 输出结果：PR1_1_fastqc.html Filename PR1_1.fq File type Conventional base calls Encoding Illumina 1.5 Total Sequences 105300 #总序列数 Sequences flagged as poor quality 0 Sequence length 90 #序列长度 %GC 52 #GC碱基含量
质量评估报告，使用浏览器打开：
3. 过滤低质量序列
利用Trimmomatic软件除去序列文件中的接头（adapter），并对碱基进行合适的修改，然后对碱基进行修剪，对低质量的序列进行过滤。
time java -jar trimmomatic-0.39.jar PE #trimmomatic是依赖java环境运行 -threads 1 -phred33 PR1_1.fq PR1_2.fq -summary ../01_trimmomaticFiltering/PR1.summary -baseout ../01_trimmomaticFiltering/PR1 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36
运行程序对序列文件中低质量的序列进行过滤，将输出结果存储到01_trimmomaticFiltering，每一个样本序列会生成如下输出文件，summary文件包含过滤的总结信息。
-rw-r--r-- 1 19M Nov 28 14:29 PR1_1P -rw-r--r-- 1 0 Nov 28 14:29 PR1_1U -rw-r--r-- 1 19M Nov 28 14:29 PR1_2P -rw-r--r-- 1 0 Nov 28 14:29 PR1_2U -rw-r--r-- 1 282 Nov 28 14:29 PR1.summary
打开其中一个PR1.summary文件进行查看，其中Input Read Pairs表示过滤之前的reads数，Both Surviving Reads表示过滤之后的reads数。
$ cat PR1.summary #查看总结文件 Input Read Pairs: 102300 Both Surviving Reads: 102300 Both Surviving Read Percent: 100.00 Forward Only Surviving Reads: 0 Forward Only Surviving Read Percent: 0.00 Reverse Only Surviving Reads: 0 Reverse Only Surviving Read Percent: 0.00 Dropped Reads: 0 Dropped Read Percent: 0.00
4. 比对到参考基因组
利用hisat2软件，将fasta序列比对到参考基因组。首先需要构建索引文件，下载或者拷贝参考基因组序列文件Ref和基因组注释文件，通过hisat2-build命令构建索引文件。
cd xx/RNA-Seq_Practice cp -r xx/Ref . cd Ref gunzip -c chr.fa.gz > chr.fa #解压参考基因组 time hisat2-build -p 1 chr.fa Chr #建立索引文件
完成上述步骤后，将过滤得到的reads比对到参考基因组上。输入文件为两个fasta序列文件，将运行过程中的输出提示重定向到log文件。
- -S设置输出文件为sam
- -x设置参考基因组
- -p设置运行的线程数量
time hisat2 -p 1 #线程数为1 -x Ref/Chr #参考基因组文件路径 -1 01_trimmomaticFiltering/PR1_1P #输入文件 -2 01_trimmomaticFiltering/PR1_2P #输入文件 -S 02_hisat2Mapping/PR1.sam #输出文件sam格式 --new-summary 1>02_hisat2Mapping/PR1_hisat2Mapping.log 2>&1 #输出日志
得到的日志文件中包含比对成功的reads数量和比对率等信息：
HISAT2 summary stats: Total pairs: 102300 Aligned concordantly or discordantly 0 time: 94209 (92.09%) Aligned concordantly 1 time: 7613 (7.44%) Aligned concordantly >1 times: 293 (0.29%) Aligned discordantly 1 time: 185 (0.18%) Total unpaired reads: 188418 Aligned 0 time: 186684 (99.08%) Aligned 1 time: 1575 (0.84%) Aligned >1 times: 159 (0.08%) Overall alignment rate: 8.76% #总比对率8.76%
上述步骤完成后，会在当前目录下生成多个sam格式文件，SAM（sequence Alignment/mapping)格式是高通量测序中存放比对数据的标准格式。另外，bam是sam的二进制格式，可以减小sam文件的大小，因此利用samtools对sam进一步处理，得到bam文件，以下用PR1为例，其他样本按相同方式处理。
time samtools view -bS 02_hisat2Mapping/PR1.sam > 02_hisat2Mapping/PR1.bam time samtools sort 02_hisat2Mapping/PR1.bam > 02_hisat2Mapping/PR1.srt.bam
5. 基因表达定量分析
利用featuresCounts软件，对reads进行精确的计数，最后将所有样本的reads数合并为一个文件，将RNA-Seq_Practice_countstable文件导出，根据FPKM和TPM的计算公式定量 每个基因在每个样本中的表达量。利用网站（https://www.ncbi.nlm.nih.gov/gene）将基因ID转化为Gene description，从而了解其功能和相关信息。
time featureCounts -a RefMaize_AGPv4/Zea_mays.AGPv4.32.gtf.gz #基因组注释文件路径 -T 1 -p --countReadPairs -g gene_id -t exon -o 03featurecountsQuatification/PR1 02hisat2Mapping/PR1.srt.bam
以PR1为例，将多余的信息除去，只保留基因ID和reads数，得到count文件，然后同样步骤处理其他样本，最后将所有的reads定量数据合并为一个countstable文件。
cat PR1 | cut -f 1,7 > PR1.count paste PR1.count PR2.count PR3.count SR1.count SR2.count SR3.count > countstable #合并不同样本的定量信息到一个文件 awk '{$3="";$5="";$7="";$9="";$11="";print $0}' countstable > RNA-Seq_Practice_countstable #提取指定列到新文件
将生成的RNA-Seq_Practice_countstable保存到本地，然后计算FPKM和TPM值，在R语言中进行相关计算。
“FPKM
”（Fragments Per Kilobase of exon model per Million mapped fragments）表示每千个碱基的转录每百万映射读取的fragments，该方法是利用每个样本的总fragments数进行校正。优点是消除样本间和基因间的差异，但是该方法主要是计算的结果可能与真实表达水平有偏差。
“TPM
”(transcript per million)表示每百万转录本中来自于某基因的转录本数目，该方法先对每个基因的reads数用基因的长度进行校正，然后再利用校正后的reads数与校正后该样本所有的reads数求商。优点是消除了外显子长度造成的差异和样本间测序总reads counts不同造成的差异，缺点为该法不是采用比对到基因组上的总reads counts，所以有时候不够精确。
具体计算方法如下所述：
> df <- read.table("RNA-Seq_Practice_countstable") #读入数据文件 > rownames(df) <- df$V1 #修改数据表的行名为基因ID（第一列V1） > df <- df[,c(2,3,4,5,6,7)] #现在行名已经为基因ID，因此删除数据的第一列 > colnames(df) <- c("PR1","PR2","PR3","SR1","SR2","SR3") #修改行名 > dim(df); names(df) #查看df的数据结构和变量名称是否正确 [1] 41768 6 [1] "PR1" "PR2" "PR3" "SR1" "SR2" "SR3" > head(df[,1:6],2) #输出前两行数据进行预览 PR1 PR2 PR3 SR1 SR2 #PR和SR共六列数据 Zm00001d027 0 0 0 0 0 Zm00001d021 0 0 0 0 0 featureCounts_meta <- df[,1:6] #提取基因信息 df <- df[rowSums(df)>0,] #去除表达量均为0的行 prefix <-"maize_PR_SR" #设置输出文件前缀名 # ----- TPM计算 ------ kb <- nrow(featureCounts_meta)/ 1000 rpk <- df / kb tpm <- t(t(rpk)/colSums(rpk) * 1000000) write.table(tpm, paste0(prefix,"_tpm.xls"), quote=F, sep="\t", row.names=T, col.names=T ) # ----- FPKM计算 ------ fpkm <- t(t(rpk)/colSums(df) * 10^6) write.table(fpkm,file= paste0(prefix, "_fpkm.xls"), quote=F, sep="\t", row.names=T, col
