轉(zhuǎn)錄組分析在當下研究功能基因組領(lǐng)域十分常用棚品。相關(guān)軟件組合種類也十分豐富，本文采用了hisat2+samtools+stringtie策略從轉(zhuǎn)錄組數(shù)據(jù)中挖掘差異表達基因。在這里小編整理了一下此套組合的執(zhí)行流程，以供日后查閱；同時分享在平臺搅轿，希望能幫助到更多初學(xué)者，如有謬誤也請各路大佬批評指正富玷。博文中涉及的參考基因組示例數(shù)據(jù)已經(jīng)上傳至我的github璧坟，可通過如下鏈接訪問下載：https://github.com/Hangyuan-Cheng/hisat2-samtools-stringtie-for-RNA-seq
先從整體上看一下軟件們所執(zhí)行的功能：
hisat2：建立參考基因組索引，reads的比對
samtools：sam2bam的轉(zhuǎn)化
stringtie：估算轉(zhuǎn)錄本表達量
所用的數(shù)據(jù)結(jié)構(gòu)如下赎懦，此處用了酵母的轉(zhuǎn)錄組數(shù)據(jù)和參考基因組：
雙端測序數(shù)據(jù)已經(jīng)經(jīng)過fastqc過濾雀鹃，具體過濾流程本文沒有涉及。示例數(shù)據(jù)僅供參考励两。執(zhí)行時要關(guān)注文件格式轉(zhuǎn)化黎茎，以及各種格式下包含的生物信息。
@biocloud:~/1223/NGS2022$ tree
.
├── gene_data.csv//差異表達基因樣例当悔，相當于最終輸出的標準答案傅瞻。
├── genome
│ ├── yeast.fa//參考基因組文件
│ ├── yeast.gff//參考基因組注釋文件
│ └── yeast_transcriptome.fa // Kallisto 所需的轉(zhuǎn)錄組索引文件，實現(xiàn)基于 pseudo alignment 的轉(zhuǎn)錄本定量時才會用到盲憎，本文不涉及該文件嗅骄。
├── reads//雙端測序clean data，即已經(jīng)經(jīng)歷過reads過濾饼疙。通常下機的時候公司都會同時給出rawdata和cleandata溺森，直接用后者就好。
│ ├── s1_y_1.fq.gz
│ ├── s1_y_2.fq.gz
│ ├── s2_y_1.fq.gz
│ └── s2_y_2.fq.gz
├── script //腳本文件，需要時加入絕對路徑引用該腳本即可屏积。
│ ├── edgeR.R
│ ├── prepDE.py
│ ├── prepDE.py3
│ └── run.sh
└── src//可能用到的軟件安裝包医窿，服務(wù)器如果安裝過該軟件則無需再次安裝。
├── fastqc_v0.11.9.zip
├── hisat2-2.2.1-Linux_x86_64.zip
├── hisat2-2.2.1-OSX_x86_64.zip
├── kallisto_linux-v0.46.1.tar.gz
├── kallisto_mac-v0.46.1.tar.gz
├── samtools-1.14.tar.bz2
├── stringtie-2.2.0.Linux_x86_64.tar.gz
└── stringtie-2.2.0.OSX_x86_64.tar.gz
1.hisat2構(gòu)建參考基因組索引
//進入存儲有參考基因組yeast.fa的目錄genome
$ cd ./genome
//將參考基因組yeast.fa建立索引并重命名為genome.fa炊林。注意hisat2與-build之間不要加空格姥卢。
$ hisat2-build yeast.fa genome.fa
2.hisat2將同一個處理下的雙端測序reads比對到參考基因組
//重新新建文件夾alignment并在其中執(zhí)行比對操作
$ cd ..
$ mkdir alignment
$ cd alignment
//hisat2將同一個處理下的測序得到的雙端測序reads比對到建立好索引的參考基因組。
//hisat2 參數(shù)-1和-2后面跟著雙端測序數(shù)據(jù)（.gz格式）路徑铛铁。比對后的結(jié)果保存為s1.sam隔显。
$ hisat2 -p 6 -x ../genome/genome.fa -1 ../reads/s1_y_1.fq.gz -2 ../reads/s1_y_2.fq.gz -S ./s1.sam
3.samtools將sam文件進行二進制轉(zhuǎn)化却妨，排序以及添加索引饵逐。
//samtools view將剛剛得到的s1.sam轉(zhuǎn)化為s1.bam
$ samtools view -bS s1.sam -o s1.bam //bam是sam的二進制文件，二進制轉(zhuǎn)化后占用存儲空間更小彪标。
//samtools sort將s1.bam排序倍权，產(chǎn)生s1.sorted.bam。后者會自動加上.bam后綴捞烟，無需命令行里添加薄声。
$ samtools sort s1.bam s1.sorted
//samtools index將排序后的bam文件添加索引
$ samtools index s1.sorted.bam
4.stringtie根據(jù)gtf注釋和排序后bam比對結(jié)果估算轉(zhuǎn)錄本表達量
//輸入剛剛排序好的s1.sorted.bam文件，酵母基因組注釋文件yeast.gff题画。
//運行stringtie得到gtf文件(s1_out.gtf)和基因表達列表(s1_genes.list)默辨。
$ stringtie ./s1.sorted.bam -G ../genome/yeast.gff -e -p 2 -o ./s1_out.gtf -A ./s1_genes.list
用同樣的方法得到另一處理下雙端測序的比對結(jié)果：s2_out.gtf 和 s2_genes.list。步驟與上述完全一致苍息，此處不再逐一重復(fù)缩幸。每個處理下的一對雙端測序結(jié)果會產(chǎn)生一個s1_out.gtf文件。按照上述步驟處理完另一對雙端測序后應(yīng)該得到兩個gtf文件：s1_out.gtf 和 s2_out.gtf竞思。
5.匯總兩個處理下的轉(zhuǎn)錄組表達量
#在新建的文件夾執(zhí)行匯總表達量操作
$ cd ..
$ mkdir differential_expression
$ cd differential_expression
#手動新建并編輯一個sample_list.txt文件表谊，用以存放處理的名稱和上述得到的gtf文件路徑。
$ vim sample_list.txt
sample_list.txt文件格式如下所示：兩個處理的名稱（s1盖喷，s2）+對應(yīng)gtf文件的地址爆办。注意第二行末尾不要有換行符號（\n）！？问帷距辆！處理名稱（s1，s2）和gtf文件地址之間應(yīng)該有一個空格暮刃。
s1 /mnt/1223/NGS2021/differential_expression/s1_out.gtf
s2 /mnt/1223/NGS2021/differential_expression/s2_out.gtf
6.使用stringtie的prepDE.py3腳本生成差異表達基因列表挑格。
prepDE.py3腳本可以在github中下載（https://github.com/gpertea/stringtie）
// 注意stringtie差異基因匯總腳本有prepDE.py和prepDE.py3兩個版本
// 前者適用于python2環(huán)境，后者python3環(huán)境沾歪。使用混了會報錯漂彤。本文基于python3.9.5環(huán)境。
$ python prepDE.py3 -i sample_list.txt -g gene_count.csv -t transcript.csv
執(zhí)行完上述步驟后得到的gene_count.csv即為差異表達基因列表，可以導(dǎo)入到R語言用edgeR / DESeq2包進行差異表達基因分析挫望，有機會再整理后繼教程立润。
sam文件基礎(chǔ)
sam文件在序列分析中至關(guān)重要，無論是轉(zhuǎn)錄組分析還是基因組call SNP媳板。認識sam文件所包含的信息有助于理解數(shù)據(jù)桑腮。我們依然以本文中生成sam文件的命令行舉例。轉(zhuǎn)錄組分析核心比對步驟是使用hisat2將測序得到的reads比對到建立好索引的參考基因組：
$ hisat2 -p 6 -x ../genome/genome.fa -1 ../reads/s1_y_1.fq.gz -2 ../reads/s1_y_2.fq.gz -S ./s1.sam
- -x 參考基因組 .fa文件
- -1 雙端測序第1段 .fq.gz格式
- -2 雙端測序第2段 .fq.gz格式
- -S 輸出文件地址+名字蛉幸，輸出結(jié)果為SAM格式
雙端測序一般為read1.fq / read1.fq.gz / read1.fq.bz2格式破讨，前后兩端同時出現(xiàn)，并同時作為hisat2輸入奕纫。將測序得到的fastq文件經(jīng)過hisat2比對到參考基因組得到SAM文件提陶。SAM文件就是序列比對文件，有上下兩個部分匹层，分別包括頭部注釋部分和比對結(jié)果部分：
頭部注釋部分
//頭部注釋部分以@開頭：
@HD VN:1.0 SO:unsorted //HD行：VN的版本以及比對排序類型隙笆。此處顯示SO:unsorted表示沒有排列順序。
@SQ SN:NC_001133.9 LN:230218//SQ行：參考序列目錄升筏。SN：參考序列名字撑柔。LN：參考序列長度
@SQ SN:NC_001148.4 LN:948066
@SQ SN:NC_001224.1 LN:85779
@SQ SN:NC_001140.6 LN:562643
@SQ SN:NC_001141.2 LN:439888
@SQ SN:NC_001142.9 LN:745751
@SQ SN:NC_001143.9 LN:666816
@PG ID:hisat2 PN:hisat2 VN:2.0.5 CL:"/mnt/bai/public/bin/hisat2-2.0.5/hisat2-align-s --wrapper basic-0 -p 8 -x ../genome/genome.fa --dta -S ./s1.sam -1 /tmp/1909596.inpipe1 -2 /tmp/1909596.inpipe2"
//PG行：使用的比對程序名，此處為hisat2
比對結(jié)果部分
比對結(jié)果部分每行表示一個read和參考基因組的比對結(jié)果信息您访。前11列為主列铅忿，包含了大部分重要信息。
SRR5511068.3 99 NC_001147.6 1002516 60 75M = 1002624 184 GTACTTAACATTCTTCTAATCATGTTAAAAGGTAAAACCTGGCCCATTTTACGATCGATCTGTAAAATCTTATAC AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEEAEEEEEEEEEEEAEEEEEEEAEEEEEEEEEEA AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YS:i:0 YT:Z:CP NH:i:1
SRR5511068.3 147 NC_001147.6 1002624 60 76M = 1002516 -184 GATAAGTAAGCAATGGTGGTAATTGCAATATTTTGCATATGTGCACGAGAAGAACTATTTTGAAGTAAGATCACTG /EEEE<EAEEEEEE6E/EEEEEEEAEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEAAAAA AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:76 YS:i:0 YT:Z:CP NH:i:1
SRR5511068.9 83 NC_001146.8 299462 60 75M1S = 299435 -103 ATTGGAAAAGAAAGTCGCGGCAAAGAGAAATGCCAATAAGACCGGGAATCAAAATTCTAAAAAGAAGAGTCAGAAG <EAAA6<A<E/EEA/A</AE/EEEAA<EEAE6AEEAAEEAEEEAAAEE/EEEEEEEEEEEE//EEEEEE/EAAAAA AS:i:-1 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YS:i:0 YT:Z:CP NH:i:1
主體部分以tab分隔從左到右依次為：
- QNAME：Query Name灵汪，測序reads名稱檀训。如果是雙端測序則會出現(xiàn)二次，兩端各比對上一次识虚。
-
FLAG：數(shù)值為2的次方數(shù)或者其加和肢扯，每個2的次方代表一種情況。詳情可以參考CSDN大佬的博客：https://blog.csdn.net/genome_denovo/article/details/78712972
或者簡書大佬：
http://www.reibang.com/p/ab133ee9712c - RENAME：Reference Name担锤，雙端測序R1比對上的參考序列名蔚晨，沒比上就是*。
- POS：position肛循，雙端測序R1起始比對上的位置序號铭腕。
- MAPQ：Mapping Quality，質(zhì)量分數(shù)多糠，越高代表越準確累舷。
- CIGAR：比對結(jié)果，M代表完全匹配夹孔。
- RNEXT：雙端測序R2端比對情況被盈，比對不上用*號析孽，比對到同一段用=號。
- MPOS：雙端測序R2端比對位置只怎。
- ISIZE：文庫插入長度袜瞬，R2端位置-R1端位置+CIGAR處的值。
- SEQ：序列信息身堡。
- QUAL：reads質(zhì)量邓尤，用ASCII碼表示。
最后總結(jié)一下：
轉(zhuǎn)錄組數(shù)據(jù)分析步驟不僅僅是比對產(chǎn)生差異表達基因贴谎，后繼還會涉及差異表達基因的顯著性分析汞扎，相關(guān)性分析，GO和KEGG分析等等擅这。本文只是總結(jié)了轉(zhuǎn)錄組分析的上游步驟澈魄。而就上游步驟而言，可用的比對軟件種類也非常多蕾哟，本文只是借助hisat2+samtools+stringtie的經(jīng)典步驟來學(xué)習(xí)轉(zhuǎn)錄組上游分析一忱。其他比對軟件日后如有涉及再行整理莲蜘。
