一、測(cè)序數(shù)據(jù)的準(zhǔn)備
新建工作目錄，RNASeq-analysis
mkdir?RNASeq-analysis
cd?RNASeq-analysis/
新建data文件夾；
mkdir data
1.1 公司返回的測(cè)序下機(jī)數(shù)據(jù)
例如：兩個(gè)樣本，三個(gè)生物學(xué)重復(fù)；
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數(shù)據(jù)庫下載數(shù)據(jù)
第一步，獲取SRA編號(hào)Accession: PRJNA******/SRP******；
第二步，下載安裝sra-tookit
conda install sra-tools
第三步，下載數(shù)據(jù)
prefetch SRR******
第四步，對(duì)數(shù)據(jù)進(jìn)行拆分
使用 SRA-toolkit 軟件包里的 fastq-dump 對(duì)數(shù)據(jù)進(jìn)行拆分，并轉(zhuǎn)換為 fastq 格式。
nohup fastq-dump --split-3 SRR****** 1>SRR******.log 2>SRR******.err &
將測(cè)序數(shù)據(jù)全部置于data文件夾，
mv Sample* data/
二、對(duì)數(shù)據(jù)進(jìn)行質(zhì)控和過濾
進(jìn)行質(zhì)控和過濾的軟件有很多，這里推薦使用國(guó)內(nèi)OpenGene 最新開發(fā)的 fastp，速度非?？欤乙徊酵瓿少|(zhì)控和過濾的流程，相當(dāng)于之前經(jīng)常用的fasqc質(zhì)控+Trimmomatic過濾，可以自動(dòng)檢測(cè)adapter 序列，無需手動(dòng)輸入，非常方便。其質(zhì)控和過濾效果也非常不錯(cuò)，可自行查閱相關(guān)資料。
fastp的運(yùn)行模式分為單端測(cè)序和雙端測(cè)序（當(dāng)然目前主流都是雙端），支持壓縮文件（file.fq.gz）及fastq輸入：
進(jìn)入工作目錄：
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
當(dāng)然，針對(duì)多個(gè)樣本，我們也可以寫個(gè)循環(huán)批量處理：
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
三、參考基因組準(zhǔn)備
主要準(zhǔn)備兩個(gè)文件，一個(gè)是基因組文件Species.genome.fasta; 一個(gè)是基因結(jié)構(gòu)注釋文件Species.genome.gff3;
新建references文件夾：
mkdir references
將參考文件置于references文件夾；
mv?Species.genome.fasta?Species.genome.gff3 ~/RNASeq-analysis/references/
將gff3文件轉(zhuǎn)換為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文件格式的區(qū)別參考：http://www.itdecent.cn/p/a27be34d335d
四、比對(duì)到參考基因組
RNASeq reads 由于不包含內(nèi)含子，所以來自外顯子邊界處的 reads被重新 mapping 回基因組
時(shí)，會(huì)被中間的內(nèi)含子分開，這種情況叫做 splice-alignment。比對(duì)的軟件有很多，在這里推薦使用Hisat2或STAR。但STAR需要更大的內(nèi)存，運(yùn)行時(shí)間也更長(zhǎng)，但準(zhǔn)確性卻相差不大。在這里使用Hisat2，使用conda安裝Hisat2。
conda install hisat2
第一步，為參考基因組構(gòu)建索引；
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
運(yùn)行hisat2-build.sh
nohup bash hisat2-build.sh &
第二步，使用Hisat2進(jìn)行比對(duì)；
新建hisat2文件夾：
mkdir hisat2
cd ~/RNASeq-analysis/hisat2
編寫shell腳本：
vim hisat2-run.sh
genome_prefix=~/RNASeq-analysis/references/Species
# 寫一個(gè)for循環(huán)進(jìn)行批量運(yùn)算
#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
運(yùn)行hisat2-run.sh
nohup bash hisat2-run.sh &
第三步，比對(duì)結(jié)果壓縮、排序及構(gòu)建索引
新建sorted目錄
mkdir sorted
進(jìn)入工作目錄
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
運(yùn)行sorted.sh
nohup bash?sorted.sh &
其中bam和sam文件均可使用IGV查看
五、計(jì)算表達(dá)量
第一步，表達(dá)定量和標(biāo)準(zhǔn)化
使用R版本Subreads軟件包中featureCounts計(jì)算表達(dá)量，其中Subread 的輸入是 bam 文件和 gtf 文
件，輸出文件描述了落在每一個(gè)基因上的 reads 數(shù)目（read counts）。read counts還需要 TPM 和
TMM標(biāo)準(zhǔn)化，這個(gè)標(biāo)準(zhǔn)化的步驟可以使用edgeR包里的fpkm函數(shù)。
新建featurecounts文件夾
mkdir?featurecounts
cd?~/RNASeq-analysis/featurecounts
編寫R腳本
vim featurecounts.R
關(guān)于R包安裝：
if (!requireNamespace("BiocManager", quietly = TRUE))
? ? install.packages("BiocManager")
BiocManager::install("limma")
BiocManager::install("Rsubread")
BiocManager::install("edgeR")
install.packages("reshape2")
我們寫一個(gè)循環(huán)運(yùn)行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數(shù)；第四列為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)
運(yùn)行R腳本
Rscript matrix_count.R
生成gene_count.csv
第1列為基因名；第2-4列為Sample1三個(gè)生物學(xué)重復(fù)對(duì)應(yīng)的count數(shù)；第5-7列為Sample2三個(gè)生物學(xué)重復(fù)對(duì)應(yīng)的count數(shù)；
編寫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)
運(yùn)行R腳本
Rscript matrix_tpm.R
生成gene_tpm.csv
第1列為基因名；第2-4列為Sample1三個(gè)生物學(xué)重復(fù)對(duì)應(yīng)的TPM值；第5-7列為Sample2三個(gè)生物學(xué)重復(fù)對(duì)應(yīng)的TPM值；
有參轉(zhuǎn)錄組的上游分析今天就寫到這里，接下來便是差異表達(dá)、后續(xù)個(gè)性化分析及可視化作圖了。
參考文獻(xiàn)：
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
