SRA—>FASTQ—>BAM—>COUNTS 這幾個(gè)步驟而已茄猫，中間穿插一些質(zhì)控的手段运褪，每個(gè)步驟選擇好合適的軟件即可洋满÷瘢可以參考：一個(gè)植物轉(zhuǎn)錄組項(xiàng)目的實(shí)戰(zhàn) http://www.bio-info-trainee.com/2809.html
這是RNA-Seq 上游分析的大致流程，比對+定量红符。當(dāng)然實(shí)驗(yàn)?zāi)康娜糁恍枰恳阎蚯啾部梢赃x擇free-alignment 的流程工具如kallisto/Salmon/Sailfish，其優(yōu)點(diǎn)是可用于RNA-seq的基因表達(dá)的快速定量预侯，但是對于小RNA和表達(dá)量低的基因分析效果并不好（2018年剛發(fā)表的一篇文章對free-alignment 的工具進(jìn)行了質(zhì)量評估致开，doi: https://doi.org/10.1101/246967）∥冢基于比對的流程双戳，比對工具也有很多選擇，如Hisat糜芳，STAR,Tophat(hista可以替代tophat),bowtie等, 還有據(jù)說速度超快的Subread飒货。根據(jù)比對工具的選擇魄衅，定量軟件也可以配套選擇，如
STAR/bowtie—RSEM塘辅；
Hisat —featureCounts/HTseq晃虫；
Subread—featureCounts。
由于后續(xù)要鑒定新的轉(zhuǎn)錄本扣墩，因此需要對轉(zhuǎn)錄本組裝哲银，組裝可以選擇cufflinks,或者stringtie。文章是用的Tophat+cufflinks呻惕。我的后續(xù)分析打算用Hisat2比對荆责，stringtie組裝，featureCounts定量亚脆，同時(shí)打算嘗試下Subread+featureCounts的流程做院。
Hisat2+stringtie+featureCounts；
Subread+featureCounts濒持。
SRA—>FASTQ
sra格式的數(shù)據(jù)需要先用fastq-dump轉(zhuǎn)換, --split-3 表示雙端測序山憨，--gzip將生成的fastq文件壓縮。
cd /pnas/fangxd_group/renyx/macaque/01rawdata
for ((i=230;i<=237;i++)) ;do /software/biosoft/software/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump --split-3 --gzip SRR4042$i.sra;done
質(zhì)控
mkdir QC
echo "fastqc started at $(date)" #輸出運(yùn)行時(shí)間
fastqc -o QC *gz
multiqc *fastqc.zip --ignore *.html
echo "fastqc finished at $(date)"
質(zhì)控檢測采用fastqc和multiqc聯(lián)合使用弥喉， multiqc有利于多個(gè)樣本同時(shí)展示質(zhì)量檢測結(jié)果郁竟。FastQC的檢測報(bào)告左側(cè)會(huì)給出測序質(zhì)量的一個(gè)summary,通常并不是所有的檢測項(xiàng)都會(huì)是綠色通過，會(huì)有一些警告和fail的項(xiàng)目由境。主要可以從Per base sequence quality 看一下測序堿基質(zhì)量棚亩，Per sequence GC content 看一下GC含量，如果實(shí)際的GC含量（紅線）出現(xiàn)雙峰虏杰，且導(dǎo)致后期的序列比對很低時(shí)讥蟆，需要引起注意，有可能是存在樣品污染纺阔；再者就是看一下Adapter是否去除及去除的是否干凈瘸彤。這里的Adapter雖然顯示通過，但是去除的并不干凈笛钝，所以后續(xù)還會(huì)進(jìn)一步去除Adapter质况。
去除接頭和低質(zhì)量值
Trimmomatic 支持多線程，處理數(shù)據(jù)速度快玻靡，主要用來去除 Illumina 平臺(tái)的 fastq 序列中的接頭结榄，并根據(jù)堿基質(zhì)量值對 fastq 進(jìn)行修剪。軟件有兩種過濾模式囤捻，分別對應(yīng) SE 和 PE 測序數(shù)據(jù)臼朗，同時(shí)支持 gzip 和 bzip2 壓縮文件。另外也支持 phred-33 和 phred-64 格式互相轉(zhuǎn)化。
Trimmomatic 過濾步驟
Trimmomatic 過濾數(shù)據(jù)的步驟與命令行中過濾參數(shù)的順序有關(guān)视哑，通常的過濾步驟如下：
ILLUMINACLIP: 過濾 reads 中的 Illumina 測序接頭和引物序列绣否，并決定是否去除反向互補(bǔ)的 R1/R2 中的 R2, 該引物序列可以在Trimmomatic軟件的安裝目錄下找到，雙端通常選擇TruSeq3-PE-2挡毅。
SLIDINGWINDOW: 從 reads 的 5’ 端開始蒜撮，進(jìn)行滑窗質(zhì)量過濾，切掉堿基質(zhì)量平均值低于閾值的滑窗慷嗜。
MAXINFO: 一個(gè)自動(dòng)調(diào)整的過濾選項(xiàng)，在保證 reads 長度的情況下盡量降低測序錯(cuò)誤率丹壕，最大化 reads 的使用價(jià)值庆械。
LEADING: 從 reads 的開頭切除質(zhì)量值低于閾值的堿基。
TRAILING: 從 reads 的末尾開始切除質(zhì)量值低于閾值的堿基菌赖。
CROP: 從 reads 的末尾切掉部分堿基使得 reads 達(dá)到指定長度缭乘。
HEADCROP: 從 reads 的開頭切掉指定數(shù)量的堿基。
MINLEN: 如果經(jīng)過剪切后 reads 的長度低于閾值則丟棄這條 reads琉用。
AVGQUAL: 如果 reads 的平均堿基質(zhì)量值低于閾值則丟棄這條 reads堕绩。
TOPHRED33: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-33。
TOPHRED64: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-64邑时。
參考：NGS 數(shù)據(jù)過濾之 Trimmomatic 詳細(xì)說明
echo "trimmomatic cut adapters started at $(date)"
for i in {230..237}
do
java -jar /software/biosoft/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 4 SRR4042$i\_1.fastq.gz SRR4042$i\_2.fastq.gz \
SRR4042$i\_paired_clean_R1.fastq.gz \
SRR4042$i\_unpair_clean_R1.fastq.gz \
SRR4042$i\_paired_clean_R2.fastq.gz \
SRR4042$i\_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/software/biosoft/software/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
done
mv *_paired_clean* ../02clean_data/
echo "trimmomatic cut adapters finished at $(date)"
質(zhì)控前后的結(jié)果
Hisat2 比對+featureCounts定量
Hisat2 比對奴紧，首先用hisat2-build 構(gòu)建index，然后比對晶丘，輸出sam格式黍氮，再用samtools將sam轉(zhuǎn)為bam格式，并排序構(gòu)建索引浅浮。
featureCounts 目前已經(jīng)整合在subread沫浆，subread是一個(gè)綜合的軟件，還包括比對的工具和對應(yīng)的R包滚秩， featrueCounts的定量效果和HTSeq差不多专执，但是featureCounts的優(yōu)點(diǎn)非常快郁油！
featureCounts下載
二進(jìn)制版本的軟件解壓即可使用
#下載subread
wget https://sourceforge.net/projects/subread/files/subread-1.5.3/subread-1.5.3-Linux-x86_64.tar.gz
tar zxvf subread-1.5.3-Linux-x86_64.tar.gz
常用參數(shù)說明：
-p 針對paired-end數(shù)據(jù)
-T 多線程數(shù)
-t 默認(rèn)將exon作為一個(gè)feature
-g 默認(rèn)是gene_id本股，從注釋文件中提取Meta-features信息用于read count
-a 基因組注釋文件，默認(rèn)是gtf
-o 輸出文件
（PS：存儲(chǔ)不夠了桐腌，后續(xù)選擇兩組數(shù)據(jù)做流程分析痊末。）
echo "hisat2 started at $(date)"
#make index
cd /pnas/fangxd_group/renyx/macaque/00ref
hisat2-build -p 8 Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa maca_index
#alignment
for i in {230..231}
do
hisat2 -t -x 00ref/maca_index -1 02clean_data/SRR4042$i\_paired_clean_R1.fastq.gz \
-2 02clean_data/SRR4042$i\_paired_clean_R2.fastq.gz \
-S 03align_out/SRR4042$i\.sam
done
#covert to bam, and sort bam
for i in {230..231}
do
samtools view -Sb SRR4042$i\.sam > SRR4042$i\.bam
samtools sort SRR4042$i\.bam -o SRR4042$i\_sorted.bam
samtools index SRR4042$i\_sorted.bam
rm *sam
done
echo "hisat2 finished at $(date)"
#featureCounts定量
cd /pnas/fangxd_group/renyx/macaque/
featureCounts=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/featureCounts
$featureCounts -p -T 6 -t exon -g gene_id -a 00ref/Macaca_mulatta.Mmul_8.0.1.91.gtf \
-o 05featurecount_out/counts.txt 03align_out/hisat2_mapping/*sorted.bam
Hisat2的比對結(jié)果：
SRR4042230和SRR4042231的比對率分別為88.02%和88.81%，總比對時(shí)間將近5小時(shí)哩掺。
featureCounts結(jié)果產(chǎn)生兩個(gè)文件：
hisat_counts.txt.summary包含一些基本的統(tǒng)計(jì)信息凿叠。
hisat_counts.txt
包含8列信息，分別是Geneid；Chr 染色體號(hào)盒件，這里會(huì)顯示多個(gè)數(shù)目蹬碧；Start; End; Strand; Length; 第7列和第8列顯示的是Counts數(shù)。
subread 比對+定量
#make subread index
echo "subread started at $(date)"
buildindex=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subread-buildindex
cd /pnas/fangxd_group/renyx/macaque/00ref
$buildindex -o maca Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa
#alignment
subjunc=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subjunc
subjunc_maca_index=/pnas/fangxd_group/renyx/macaque/00ref/maca
cd /pnas/fangxd_group/renyx/macaque/
for i in {230..231}
do
$subjunc -T 5 -i $subjunc_maca_index -r 02clean_data/SRR4042$i\_paired_clean_R1.fastq.gz -R 02clean_data/SRR4042$i\_paired_clean_R2.fastq.gz -o 03align_out/subread_mapping/SRR4042$
done
#counts expre
featureCounts=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/featureCounts
$featureCounts -p -T 6 -t exon -g gene_id -a 00ref/Macaca_mulatta.Mmul_8.0.1.91.gtf \
-o 05featurecount_out/subread_counts.txt 03align_out/subread_mapping/*_subjunc.bam
echo "subread finished at $(date)"
可以看出subjunc比對不到2小時(shí)炒刁，比hisat2快3個(gè)多小時(shí)恩沽。
提取counts
根據(jù)第1列是Geneid，第7,8列是counts數(shù)翔始，用awk提取出geneID和counts罗心。
awk -F '\t' '{print $1,$7,$8}' OFS='\t' hisat_counts.txt >hisat_matrix.out
awk -F '\t' '{print $1,$7,$8}' OFS='\t' subread_counts.txt >subread_matrix.out
參考資料：
一個(gè)植物轉(zhuǎn)錄組項(xiàng)目的實(shí)戰(zhàn)
史上最快的轉(zhuǎn)錄組流程-subread
轉(zhuǎn)錄組定量-FEATURECOUNTS
