這是RNA-Seq 上游分析的大致流程，比對+定量。當然實驗目的若只需要定量已知基因，也可以選擇free-alignment 的流程工具如kallisto/Salmon/Sailfish，其優點是可用於RNA-seq的基因表達的快速定量，但是對於小RNA和表達量低的基因分析效果並不好（2018年剛發表的一篇文章對free-alignment 的工具進行了質量評估，doi: https://doi.org/10.1101/246967）。基於比對的流程，比對工具也有很多選擇，如Hisat，STAR,Tophat(hista可以替代tophat),bowtie等, 還有據說速度超快的Subread。根據比對工具的選擇，定量軟體也可以配套選擇，如 STAR/bowtie—RSEM； Hisat —featureCounts/HTseq； Subread—featureCounts。 由於後續要鑑定新的轉錄本，因此需要對轉錄本組裝，組裝可以選擇cufflinks,或者stringtie。文章是用的Tophat+cufflinks。我的後續分析打算用Hisat2比對，stringtie組裝，featureCounts定量，同時打算嘗試下Subread+featureCounts的流程。 Hisat2+stringtie+featureCounts；

Subread+featureCounts。

SRA—>FASTQ

sra格式的資料需要先用fastq-dump轉換, --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

質控

mkdir QC
echo "fastqc started at $(date)"     #輸出執行時間
fastqc -o QC *gz
multiqc *fastqc.zip --ignore *.html
echo "fastqc finished at $(date)"
 

質控檢測採用fastqc和multiqc聯合使用， multiqc有利於多個樣本同時展示質量檢測結果。FastQC的檢測報告左側會給出測序質量的一個summary,通常並不是所有的檢測項都會是綠色通過，會有一些警告和fail的專案。主要可以從Per base sequence quality 看一下測序鹼基質量，Per sequence GC content 看一下GC含量，如果實際的GC含量（紅線）出現雙峰，且導致後期的序列比對很低時，需要引起注意，有可能是存在樣品汙染；再者就是看一下Adapter是否去除及去除的是否乾淨。這裡的Adapter雖然顯示通過，但是去除的並不乾淨，所以後續還會進一步去除Adapter。

Summary QC

Adapter Content

去除接頭和低質量值

Trimmomatic 支援多執行緒，處理資料速度快，主要用來去除 Illumina 平臺的 fastq 序列中的接頭，並根據鹼基質量值對 fastq 進行修剪。軟體有兩種過濾模式，分別對應 SE 和 PE 測序資料，同時支援 gzip 和 bzip2 壓縮檔案。另外也支援 phred-33 和 phred-64 格式互相轉化。

Trimmomatic 過濾步驟

Trimmomatic 過濾資料的步驟與命令列中過濾引數的順序有關，通常的過濾步驟如下： ILLUMINACLIP: 過濾 reads 中的 Illumina 測序接頭和引物序列，並決定是否去除反向互補的 R1/R2 中的 R2, 該引物序列可以在Trimmomatic軟體的安裝目錄下找到，雙端通常選擇TruSeq3-PE-2。 SLIDINGWINDOW: 從 reads 的 5’ 端開始，進行滑窗質量過濾，切掉鹼基質量平均值低於閾值的滑窗。 MAXINFO: 一個自動調整的過濾選項，在保證 reads 長度的情況下儘量降低測序錯誤率，最大化 reads 的使用價值。 LEADING: 從 reads 的開頭切除質量值低於閾值的鹼基。 TRAILING: 從 reads 的末尾開始切除質量值低於閾值的鹼基。 CROP: 從 reads 的末尾切掉部分鹼基使得 reads 達到指定長度。 HEADCROP: 從 reads 的開頭切掉指定數量的鹼基。 MINLEN: 如果經過剪下後 reads 的長度低於閾值則丟棄這條 reads。 AVGQUAL: 如果 reads 的平均鹼基質量值低於閾值則丟棄這條 reads。 TOPHRED33: 將 reads 的鹼基質量值體系轉為 phred-33。 TOPHRED64: 將 reads 的鹼基質量值體系轉為 phred-64。 參考：NGS 資料過濾之 Trimmomatic 詳細說明

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)"

質控前後的結果

General Statistics

Hisat2 比對+featureCounts定量

Hisat2 比對，首先用hisat2-build 構建index，然後比對，輸出sam格式，再用samtools將sam轉為bam格式，並排序構建索引。 featureCounts 目前已經整合在subread，subread是一個綜合的軟體，還包括比對的工具和對應的R包， featrueCounts的定量效果和HTSeq差不多，但是featureCounts的優點非常快！

featureCounts下載

二進位制版本的軟體解壓即可使用

#下載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

常用引數說明：

-p 針對paired-end資料 T 多執行緒數 t 預設將exon作為一個feature g 預設是gene_id，從註釋檔案中提取Meta-features資訊用於read count a 基因組註釋檔案，預設是gtf o 輸出檔案

（PS：儲存不夠了，後續選擇兩組資料做流程分析。）

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的比對結果：

SRR4042230和SRR4042231的比對率分別為88.02%和88.81%，總比對時間將近5小時。

Hisat2 mapping

featureCounts結果產生兩個檔案：

hisat_counts.txt.summary包含一些基本的統計資訊。

hisat_counts.txt.summary

hisat_counts.txt包含8列資訊，分別是Geneid；Chr 染色體號，這裡會顯示多個數目；Start; End; Strand; Length; 第7列和第8列顯示的是Counts數。

hisat_counts.txt

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)"

image.png

image.png

可以看出subjunc比對不到2小時，比hisat2快3個多小時。

提取counts

根據第1列是Geneid，第7,8列是counts數，用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

參考資料：

一個植物轉錄組專案的實戰 史上最快的轉錄組流程-subread 轉錄組定量-FEATURECOUNTS

後續分析，請大家持續關注~