最近在学习转录组数据的整套分析流程，顺便记录一下。
转录组是在特定时空条件下细胞中基因转录表达产物，广义的转录组包括信使RNA，核糖体RNA，转运RNA及非编码RNA，狭义上是指所有mRNA的集合，转录组分析能够获得不同基因的表达情况。
步骤
1.数据来源
这里使用的是茶树不同组织的样本，共6个组织，每个组织三个生物学重复，共18个样本。利用illumina平台进行转录组测序，得到双端测序数据。数据原始格式为.fq,共有12条测序数据文件（每个样本产生两条）
BUD1_1.fq.gz BUD1_2.fq.gz BUD2_1.fq.gz BUD2_2.fq.gz BUD3_1.fq.gz BUD3_2.fq.gz
FLOWER1_1.fq.gz FLOWER1_2.fq.gz FLOWER2_1.fq.gz FLOWER2_2.fq.gz FLOWER3_1.fq.gz FLOWER3_2.fq.gz
MATURE1_1.fq.gz MATURE1_2.fq.gz MATURE2_1.fq.gz MATURE2_2.fq.gz MATURE3_1.fq.gz MATURE3_2.fq.gz
ROOT1_1.fq.gz ROOT1_2.fq.gz ROOT2_1.fq.gz ROOT2_2.fq.gz ROOT3_1.fq.gz ROOT3_2.fq.gz
STEM1_1.fq.gz STEM1_2.fq.gz STEM2_1.fq.gz STEM2_2.fq.gz STEM3_1.fq.gz STEM3_2.fq.gz
TENDER1_1.fq.gz TENDER1_2.fq.gz TENDER2_1.fq.gz TENDER2_2.fq.gz TENDER3_1.fq.gz TENDER3_2.fq.gz
创建工作目录
mkdir RNA-Seq_Practice
cd RNA-Seq_Practice
mkdir 00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_stringtie 04_Ballgown 05_TransDecoder
2. 测序数据质量评估
需要先安装fastQC软件，利用fastQC对获得的fastq序列文件进行质量分析，生成html格式的结果报告，包含各项指标。fastQC是一个java软件，需要自行配制JAVA环境。这里使用java -version查看java安装情况。若没有安装，请参考https://zhuanlan.zhihu.com/p/28924831
#fastQC的安装
#通过conda安装
conda install fastqc
#通过apt安装
sudo apt install fastqc
#通过源码安装
wget -c <https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip>
unzip fastqc_v0.12.1.zip
sudo ln -s /path/to/FastQC/fastqc /usr/local/bin/fastqc #创建fastqc的软链接
下载好软件后，进行质控分析，这里以BUD1为例。
fastqc *.fq.gz #切换到rawdata所在文件夹内，批量对fq后缀文件运行fastqc程序
输出结果：BUD1_1_fastqc.html
Filename BUD1_1.fq.gz
File type Conventional base calls
Encoding Sanger/Illumina 1.9
Total Sequences 24141589 #总序列数
Sequences flagged as poor quality 0
Sequence length 150 #测序长度
%GC 45 #GC碱基含量
质量评估报告，可使用浏览器打开查看。
3. 过滤低质量序列
这里使用Trimmomatic软件进行低质量序列的过滤，主要包括去除序列文件中的接头（adapter），并对碱基进行合适的修改和修剪。
安装
这里介绍两种方法，分别是conda安装和网站直接下载
#conda安装
pip install --upgrade -i https://pypi.doubanio.com/simple/argparse
conda install -c bioconda trimmomatic
#帮助文档
trimmomatic -h
#网站下载
#下载
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
#解压
unzip Trimmomatic-0.39.zip
下载后运行
time java -jar trimmomatic-0.39.jar PE #trimmomatic是依赖java环境运行
-threads 1
-phred33 BUD1_1.fq BUD1_2.fq
-summary ../01_trimmomaticFiltering/BUD1.summary
-baseout ../01_trimmomaticFiltering/BUD1 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36
具体参数见https://www.cnblogs.com/Sunny-King/archive/2022/06/16/Bioinformatics-Trimmomatic.html
输出结果存储到01_trimmomaticFiltering，每一个样本序列会生成4个输出文件，summary文件包含过滤的总结信息。
打开其中一个summary文件，可以看到具体信息.其中Input Read Pairs表示过滤之前的reads数，Both Surviving Reads表示过滤之后的reads数。
$ cat BUD1.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序列比对到参考基因组。
Hisat2安装
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip
echo 'export PATH=/home/gfq/hiasat2/hisat2-2.1.0/bin:$PATH' >> ~/.bashrc#这里若不配置环境变量则需要加上绝对路径运行
source ~/.bashrc
首先需要构建索引文件，下载参考基因组序列文件fasta和基因组注释文件，通过hisat2-build命令构建索引文件。
cd xx/Tieguanyin_Ref
gunzip -c chr.fa.gz > Tieguanyin.fa #解压参考基因组序列文件
time hisat2-build -p 1 Tieguanyin.fa Chr #建立索引文件
将过滤得到的reads比对到参考基因组上
time hisat2 -p 1 #线程数为1
-x Tieguanyin_Ref/Chr #参考基因组文件路径
-1 01_trimmomaticFiltering/BUD1_1P #输入文件
-2 01_trimmomaticFiltering/BUD1_2P #输入文件
-S 02_hisat2Mapping/BUD1.sam #输出文件sam格式
--new-summary 1>02_hisat2Mapping/BUD1_hisat2Mapping.log 2>&1
#将运行过程中的输出提示重定向到log文件，输出日志。
得到的日志文件中包含比对成功的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文件。
安装
#conda安装
conda install samtools
#普通下载安装
wget -c https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar jxvf samtools-1.9.tar.bz2#解压
cd samtools-1.9/
./configure --prefix=/home/gfq/samtools/samtools-1.9
make
make install
echo 'export PATH="/home/gfq/samtools/samtools-1.9/bin:$PATH" ' >>~/.bashrc
source ~/.bashrc
SAM文件转Bam文件，并对bam 文件中的内容进行排序。以下用BUD1为例，其他样本按相同方式处理。
time samtools view -bhS
02_hisat2Mapping/BUD1.sam > 02_hisat2Mapping/BUD1.bam
time samtools sort
02_hisat2Mapping/BUD1.bam > 02_hisat2Mapping/BUD1.srt.bam
5.利用StringTie进行转录本组装、量化基因表达
安装
wget http://ccb.jhu.edu/software/stringtie/dl/stringtie-2.0.4.Linux_x86_64.tar.gz
tar -zvxf stringtie-2.0.4.Linux_x86_64.tar.gz
echo 'export PATH=/home/gfq/stringtie/stringtie-2.0.4.Linux_x86_64:$PATH' >> ~/.bashrc
source ~/.bashrc
1.样本组装
比对上的reads将会被呈递给StringTie进行转录本组装，StringTie单独的对每个样本进行组装，在组装的过程中估算每个基因及isoform的表达水平。
stringtie -p 4 -G /home/gfq/Tieguanyin_Ref/Tieguanyin.genome.gtf -o BUD1.gtf -l BUD1 BUD1.srt.bam
#-p 线程（CPU)数 (default: 1)
#-G 参考序列的基因注释文件 (GTF/GFF3)
#-l 输出转录本的名称前缀（默认为MSTRG）
2.转录本组装
所有转录本都被呈递给StringTie的merge函数进行merge，这一步是必须的，因为有些样本的转录本可能仅仅被部分reads覆盖，无法被StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本。
vi mergelist.txt #需要包含之前output.gtf文件的路径，将所有上一步输出的名称填入文件。这里只列举BUD样本，需要把所有样本都写入。
BUD1.gtf
BUD2.gtf
BUD3.gtf
3.转录本合并
stringtie --merge -p 4 -G /home/gfq/Tieguanyin_Ref/Tieguanyin.genome.gtf -o stringtie_merged.gtf mergelist.txt
#-p 线程（CPU)数 (default: 1)
#-G <guide_gff> 参考转录本的注释信息 (GTF/GFF3)
4.检测相对于注释基因组转录本的组装情况
使用gffCompare实用程序来确定有多少组合的转录本完全或部分匹配带注释的基因，并计算出有多少是完全新颖的
安装
wget http://ccb.jhu.edu/software/stringtie/dl/gffcompare-0.11.5.Linux_x86_64.tar.gz
tar -zxvf gffcompare-0.11.5.tar.gz
echo 'export PATH=/home/gfq/gffcompare/gffcompare-0.11.5.Linux_x86_64:$PATH' >> ~/.bashrc
source ~/.bashrc
运行
gffcompare -r /home/gfq/Tieguanyin_Ref/Tieguanyin.genome.gtf -G -o merged stringtie_merged.gtf
#-r 参考转录本的注释信息
#-G 比对所有转录本
#-o 指定要用于gffcompare将创建的输出文件的前缀
结果会生成6个文件
gffcompare结果
gffcmp.annotated.gtf：这里面向每个转录本添加一个"类代码"和来自参考注释文件的转录本的名称，这使用户能够快速检查预测的转录本与参考基因组的关系。
gffcmp.stats：包含不同基因特征的灵敏度和精度
灵敏度：参考基因组中正确重建的的基因比例
精度：显示与参考基因组重叠的gene的比例
5.重新组装转录本并估算基因表达丰度
mkdir ballgown
cd 03_stringtie/
stringtie –e –B -p 4 -G stringtie_merged.gtf -o /home/gfq/RNA-Seq/04_ballgown/BUD1/BUD1.gtf BUD1.srt.bam
#-e 只对参考转录本进行丰度评估 (requires -G)
#-G 参考序列的基因注释文件 (GTF/GFF3)
#-B 在输出的GFT同目录下输出Ballgown table 文件
6.stringtie输出的表达量结果转换为表格
#下载脚本
wget http://ccb.jhu.edu/software/stringtie/dl/prepDE.py
#转换格式
python2 ./prepDE.py - ballgown
这里会生成gene_count_matrix.csv、transcript_count_matrix.csv文件。若想得到TPM/FPKM表格，也可下载对应的.py文件进行转化。详见https://www.jianshu.com/p/53a13af6f51f
6.TransDecoder预测CDS
TransDecoder能够从转录本序列中鉴定候选编码区。这些转录本序列可以来自于Trinity的从头组装或者来自于Cufflinks或者StringTie的组装结果。
安装
mkdir TransDecoder
wget https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.5.0.zip
unzip TransDecoder-v5.5.0.zip
mv TransDecoder-TransDecoder-v5.5.0 TransDecoder-v5.5.0
输入文件
transcripts.gtf: 记录预测转录本的GTF文件
genome.fasta: 参考基因组序列
从GTF文件中提取FASTA序列
cd 05_TransDecoder/
/home/gfq/TransDecoder/TransDecoder-v5.5.0/util/gtf_genome_to_cdna_fasta.pl /home/gfq/RNA-Seq/03_stringtie/transcripts.gtf /home/gfq/Tieguanyin_Ref/Tieguanyin.fasta > transcripts.fasta
将GTF文件转成GFF3格式
/home/gfq/TransDecoder/TransDecoder-v5.5.0/util/gtf_to_alignment_gff3.pl transcripts.gtf > transcripts.gff3
预测转录本中长的开放阅读框, 默认是100个氨基酸，可以用-m修改
/home/gfq/TransDecoder/TransDecoder-v5.5.0/TransDecoder.LongOrfs -t transcripts.fasta
使用DIAMOND（详情请看https://xuzhougeng.top/archives/Fast-and-sensitive-protein-alignment-using-diamond）对上一步输出的transcripts.fasta.transdecoder.pep在蛋白数据库中进行搜索，寻找同源证据支持
#安装DIAMOND
wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
mv diamond ~/bin
# 下载uniprot数据并解压缩
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip uniprot_sprot.fasta.gz
# 建立索引
diamond makedb --in uniprot_sprot.fasta --db uniprot_sprot.fasta
# BLASTP比对
diamond blastp -d uniprot_sprot.fasta -q transcripts.fasta.transdecoder_dir/longest_orfs.pep --evalue 1e-5 --max-target-seqs 1 > blastp.outfmt6
预测可能的编码区
/home/gfq/TransDecoder/TransDecoder-v5.5.0/TransDecoder.Predict \
-t transcripts.fasta \
--retain_blastp_hits blastp.outfmt6
生成基于参考基因组的编码区注释文件
/home/gfq/TransDecoder/TransDecoder-v5.5.0/util/cdna_alignment_orf_to_genome_orf.pl \
transcripts.fasta.transdecoder.gff3 \
transcripts.gff3 \
transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3
最终输出的文件有:
transcripts.fasta.transdecoder.pep: 最终预测的CDS对应的蛋白序列
transcripts.fasta.transdecoder.cds: 最终预测的CDS序列
transcripts.fasta.transdecoder.gff3: 最终ORF对应的GFF3
transcripts.fasta.transdecoder.bed: 以BED格式存放ORF位置信息
transcripts.fasta.transdecoder.genome.gff3: 基于参考基因组的GFF3文件
7.基因功能注释
基因功能注释方面主要参考：https://yanzhongsino.github.io/2021/05/17/omics_genome.annotation_function/
