RNA-Seq：从fastq到表达矩阵
一、测序数据的准备
新建工作目录，RNASeq-analysis
mkdir RNASeq-analysis
cd RNASeq-analysis/
新建data文件夹；
mkdir data
1.1 公司返回的测序下机数据
例如：两个样本，三个生物学重复；
Raw data: Sample1-1_R1.fq.gz; Sample1-1_R2.fq.gz; Sample1-2_R1.fq.gz; Sample1-2_R2.fq.gz; Sample1-3_R1.fq.gz; Sample1-3_R2.fq.gz; Sample2-1_R1.fq.gz; Sample2-1_R2.fq.gz; Sample2-2_R1.fq.gz; Sample2-2_R2.fq.gz; Sample2-3_R1.fq.gz; Sample2-3_R2.fq.gz;
1.2 从SRA数据库下载数据
第一步，获取SRA编号Accession: PRJNA******/SRP******；
第二步，下载安装sra-tookit
conda install sra-tools
第三步，下载数据
prefetch SRR******
第四步，对数据进行拆分
使用 SRA-toolkit 软件包里的 fastq-dump 对数据进行拆分，并转换为 fastq 格式。
nohup fastq-dump --split-3 SRR****** 1>SRR******.log 2>SRR******.err &
将测序数据全部置于data文件夹，
mv Sample* data/
二、对数据进行质控和过滤
进行质控和过滤的软件有很多，这里推荐使用国内OpenGene 最新开发的 fastp，速度非常快，且一步完成质控和过滤的流程，相当于之前经常用的fasqc质控+Trimmomatic过滤，可以自动检测adapter 序列，无需手动输入，非常方便。其质控和过滤效果也非常不错，可自行查阅相关资料。
fastp的运行模式分为单端测序和双端测序（当然目前主流都是双端），支持压缩文件（file.fq.gz）及fastq输入：
进入工作目录：
cd ~/RNASeq-analysis/data/
fastp (single -end, SE)
fastp -I SRR******.fastq -O SRR******_clean.fastq
fastp (paired -end, PE)
fastp -i Sample1-1_R1.fq.gz -o Sample1-1_R1.clean.fq.gz -I Sample1-1_R2.fq.gz -O Sample1-1_R2.clean.fq.gz
fastp -i Sample1-1_R1.fq.gz -I Sample1-1_R2.fq.gz -o Sample1-1_R1.clean.fq.gz -O Sample1-1_R2.clean.fq.gz
当然，针对多个样本，我们也可以写个循环批量处理：
for file in `ls *_R1.fq.gz | perl -lpe 's/_R1.fq.gz//'`; do fastp -i ${file}_R1.fq.gz -I ${file}_R2.fq.gz -o ${file}_R1.clean.fq.gz -O ${file}_R2.clean.fq.gz; done
三、参考基因组准备
主要准备两个文件，一个是基因组文件Species.genome.fasta; 一个是基因结构注释文件Species.genome.gff3;
新建references文件夹：
mkdir references
将参考文件置于references文件夹；
mv Species.genome.fasta Species.genome.gff3 ~/RNASeq-analysis/references/
将gff3文件转换为gtf文件：
安装gffread：conda install gffread
# gff to gtf
gffread Species.genome.gff3 -T -o Species.genome.gtf
# gtf to gff
gffread Species.genome.gtf -o- > Species.genome.gff3
gff3与gtf文件格式的区别参考：https://www.jianshu.com/p/a27be34d335d
四、比对到参考基因组
RNASeq reads 由于不包含内含子，所以来自外显子边界处的 reads被重新 mapping 回基因组
时，会被中间的内含子分开，这种情况叫做 splice-alignment。比对的软件有很多，在这里推荐使用Hisat2或STAR。但STAR需要更大的内存，运行时间也更长，但准确性却相差不大。在这里使用Hisat2，使用conda安装Hisat2。
conda install hisat2
第一步，为参考基因组构建索引；
cd ~/RNASeq-analysis/references/
编写shell脚本：
vim hisat2-build.sh
genome=~/RNASeq-analysis/references/Species.genome.fasta
genome_prefix=~/RNASeq-analysis/references/Species
hisat2-build $genome $genome_prefix 1>hisat2-build.log 2>hisat2-build.err
运行hisat2-build.sh
nohup bash hisat2-build.sh &
第二步，使用Hisat2进行比对；
新建hisat2文件夹：
mkdir hisat2
cd ~/RNASeq-analysis/hisat2
编写shell脚本：
vim hisat2-run.sh
genome_prefix=~/RNASeq-analysis/references/Species
# 写一个for循环进行批量运算
#for paired -end, PE
for sample in `ls ~/RNASeq-analysis/data/*_R1.clean.fq.gz | perl -lpe 's/_R1.clean.fq.gz//'`; do hisat2 --new-summary -p 2 -x $genome_prefix -1 ${sample}_R1.clean.fq.gz -2 ${sample}_R2.clean.fq.gz -S ~/RNASeq-analysis/hisat2/${sample}.sam 2> ~/RNASeq-analysis/hisat2/${sample}.err; done
#for single -end, SE
for sample in `ls ~/RNASeq-analysis/data/SRR*.clean.fq.gz | perl -lpe 's/SRR*.clean.fq.gz//'`; do hisat2 --new-summary -p 2 -x $genome_prefix -U $left -S ~/RNASeq-analysis/hisat2/${sample}.sam 2> ~/RNASeq-analysis/hisat2/${sample}.err; done
运行hisat2-run.sh
nohup bash hisat2-run.sh &
第三步，比对结果压缩、排序及构建索引
新建sorted目录
mkdir sorted
进入工作目录
cd ~/RNASeq-analysis/hisat2
编写shell脚本
vim sorted.sh
for file in *.sam; do
samtools view -b $file > ~/RNASeq-analysis/sorted/${file%.sam}.bam
samtools sort ~/RNASeq-analysis/sorted/${file%.sam}.bam > ~/RNASeq-analysis/sorted/${file%.sam}.sorted.bam
samtools index ~/RNASeq-analysis/sorted/${file%.sam}.sorted.bam
done
运行sorted.sh
nohup bash sorted.sh &
其中bam和sam文件均可使用IGV查看
五、计算表达量
第一步，表达定量和标准化
使用R版本Subreads软件包中featureCounts计算表达量，其中Subread 的输入是 bam 文件和 gtf 文
件，输出文件描述了落在每一个基因上的 reads 数目（read counts）。read counts还需要 TPM 和
TMM标准化，这个标准化的步骤可以使用edgeR包里的fpkm函数。
新建featurecounts文件夹
mkdir featurecounts
cd ~/RNASeq-analysis/featurecounts
编写R脚本
vim featurecounts.R
关于R包安装：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("limma")
BiocManager::install("Rsubread")
BiocManager::install("edgeR")
install.packages("reshape2")
我们写一个循环运行featurecounts.R得到count文件
for file in ~/RNASeq-analysis/sorted/*.sorted.bam; do Rscript featurecounts.R ~/RNASeq-analysis/sorted/$file ~/RNASeq-analysis/references/Species.genome.gtf 10 $file; done
第二步，将所有count文件合并
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count > merge.count
第一列为样品名；第二列为基因名；第三列为counts数；第四列为fpkm值；第五列为tpm值
第三步，提取counts值及tpm值并生成矩阵
awk -F"\t" '{print $1"\t"$2"\t"$3}' merge.count > gene.count
awk -F"\t" '{print $1"\t"$2"\t"$5}' merge.count > gene.tpm
编写R脚本：matrix_count.R
a=read.csv('gene.count',header=F,sep="\t")
colnames(a)=c('sample','gene','counts')
library(reshape2)
counts=dcast(a,formula=gene~sample)
write.table(counts,file="gene_count.csv",sep="\t",quote=FALSE,row.names=FALSE)
运行R脚本
Rscript matrix_count.R
生成gene_count.csv
第1列为基因名；第2-4列为Sample1三个生物学重复对应的count数；第5-7列为Sample2三个生物学重复对应的count数；
编写R脚本：matrix_tpm.R
a=read.csv('tpm.count',header=F,sep="\t")
colnames(a)=c('sample','gene','tpm')
library(reshape2)
counts=dcast(a,formula=gene~sample)
write.table(counts,file="gene_tpm.csv",sep="\t",quote=FALSE,row.names=FALSE)
运行R脚本
Rscript matrix_tpm.R
生成gene_tpm.csv
第1列为基因名；第2-4列为Sample1三个生物学重复对应的TPM值；第5-7列为Sample2三个生物学重复对应的TPM值；
有参转录组的上游分析今天就写到这里，接下来便是差异表达、后续个性化分析及可视化作图了。
参考文献：
1. Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, https://doi.org/10.1093/bioinformatics/bty560
2. Kim D, Langmead B, and Salzberg SL HISAT: a fast spliced aligner with low memory requirements, Nature methods, 2015
3. Li H.*, Handsaker B.*, Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009) The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 25, 2078-9. [PMID: 19505943]
4. Li H A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 Nov 1;27(21):2987-93. Epub 2011 Sep 8. [PMID: 21903627]
5. Liao Y, Smyth GK and Shi W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Research, 47(8):e47, 2019
6. Liao Y, Smyth GK and Shi W. featureCounts: an efficient general-purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7):923-30, 2014
7. Liao Y, Smyth GK and Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Research, 41(10):e108, 2013
8. Wickham, H. (2007). Reshaping data with the reshape package. 21(12):1–20. http://www.jstatsoft.org/v21/i12
