本節(jié)概覽：
- hisat2 + featureCounts:
獲取hisat2索引文件孽文，hisat2比對和samtools格式轉(zhuǎn)化，featureCounts計數(shù)得到counts表達矩陣- Salmon:
salmon index 用cdna.fa建立索引茵瀑，salmon quant對clean fastq文件直接進行基因定量- 獲取ensembl_id或transcript_id轉(zhuǎn)化的對應文件
承接上節(jié)RNA-seq入門實戰(zhàn)（一）
本節(jié)介紹使用hisat2或salmon這兩種方法進行轉(zhuǎn)錄組上游數(shù)據(jù)的比對和計數(shù)。
39個轉(zhuǎn)錄組分析工具，120種組合評估
(https://www.nature.com/articles/s41467-017-00050-4)
表明基于hisat2或salmon進行轉(zhuǎn)錄本定量都比較優(yōu)秀鸵荠。
一簸呈、hisat2 + featureCounts
1. 獲取hisat2比對索引文件
index官網(wǎng)下載地址Download | HISAT2 (daehwankimlab.github.io)榕订，下載并解壓所需的 mm10 或 grcm38 的index文件
mkdir -p ~/reference/index/hisat/
cd ~/reference/index/hisat/
wget https://genome-idx.s3.amazonaws.com/hisat/mm10_genome.tar.gz
tar -zxvf *tar.gz
rm *tar.gz
2. hisat2比對和samtools轉(zhuǎn)化格式
先用hisat2比對基因組得到sam文件，再用samtools sort將sam文件格式轉(zhuǎn)化與排序為bam文件（bam相當于二進制版的sam）蜕便，之后samtools index建立索引（用于后續(xù)IGV內(nèi)可視化）劫恒，最后samtools flag 統(tǒng)計文件比對情況保存在文本中。其中samtools index與samtools flag為非必須步驟轿腺，可略過两嘴。sam相當于是中間文件比較占存儲空間，可在轉(zhuǎn)化為bam后便刪除族壳。
代碼如下：
vim 3_align2sam2bam_hisat2.sh
############################
echo -e " \n \n \n 333# Align !!! hisat2 !!!\n \n \n "
date
########set#### ###set#### ###set####
index='/home/reference/index/hisat/mm10/genome'
mkdir -p ~/test/align/flag
cd ~/test/align/
pwd
cat ~/test/idname | while read id
do
echo "333# ${id} ${id} ${id} is on the hisat2 Working !!!"
################ paired ###############################
hisat2 -t -p 12 -x $index \
-1 ~/test/clean/${id}_*1*gz \
-2 ~/test/clean/${id}_*2*gz -S ${id}.sam
######################################################
################Single################################
# hisat2 -t -p 12 -x $index \
# -U ~/test/clean/${id}_trimmed.fq.gz \
# -S ./${id}.sam
########################################################
### sam2bam and remove sam ###
echo -e " ${id} sam2bam and remove sam "
samtools sort -@ 12 -o ~/test/align/${id}_sorted.bam ${id}.sam
rm ${id}.sam
done
#### samtools index and flagstat ####
echo -e " \n \n \n samtools index and flagstat \n "
cd ~/test/align/flag
pwd
#ls ~/test/align/*.bam | xargs -i samtools index -@ 12 {}
ls ~/test/align/*.bam | while read id ;\
do
samtools flagstat -@ 12 $id > $(basename $id '.bam').flagstat
done
multiqc ./
echo -e " \n \n \n 333# ALL Work Done !!!\n \n \n "
date
nohup bash 3_align2sam2bam_hisat2.sh >log_3 2>&1 &
比對結(jié)果如下：
3. featureCounts 計數(shù)得到counts表達矩陣
計數(shù)首先要獲取gtf注釋文件溶诞，注意要和hisat2的index文件的基因組版本相對應，如本次為mm10决侈，則gtf文件也必須為mm10或grcm38螺垢。
研究人和鼠推薦用gencode數(shù)據(jù)庫的文件GENCODE 喧务，比較常用的還有UCSC的refGene.gtf文件，下載地址在https://hgdownload.soe.ucsc.edu/ （若想下載其他gtf文件則將網(wǎng)址中mm10替換即可枉圃，如hg38）功茴。
featureCounts詳細使用方法見轉(zhuǎn)錄組定量可以用替換featureCounts代替HTSeq-count (qq.com)，常用參數(shù)如下：
代碼如下：
vim 4_counts_featurecounts.sh
###########################################
echo -e "\n \n \n444# Count #featureCounts !!! \n \n \n"
date
#####set####set###set
gtf='/home/reference/gtf/gencode/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf'
#gtf='/home/reference/gtf/UCSC/mm10.refGene.gtf.gz'
mkdir -p ~/test/counts
cd ~/test/counts/
pwd
######## single##########################################################################
#featureCounts -T 12 -a $gtf -o counts.txt ~/test/align/*.bam
#######paired###########################################################################
featureCounts -T 12 -p -a $gtf -o counts.txt ~/test/align/*.bam
#######################################################################################
###生成網(wǎng)頁版統(tǒng)計情況
multiqc ./
echo -e " \n \n \n ALL WORK DONE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \n "
date
nohup bash 4_counts_featurecounts.sh >log_4 2>&1 &
查看計數(shù)的統(tǒng)計情況孽亲，匹配率在71-80%左右坎穿，還可以
運行結(jié)果保存在counts文件夾下，里面的counts.txt就是我們下游分析所需要的文件啦
二返劲、Salmon——直接對基因進行定量的工具
與hisat2不同玲昧，salmon不經(jīng)過比對計數(shù)步驟而是直接對基因進行定量，如果不研究新轉(zhuǎn)錄本篮绿，用salmon方法可以更快更方便得到基因定量信息孵延。
1. 建立salmon索引
先下載參考轉(zhuǎn)錄本序列cDNA.fa文件，在ensembl官網(wǎng)選擇相應文件 Index of /pub/release-102/fasta/mus_musculus/cdna/ (ensembl.org)亲配，使用salmon index 建立索引文件尘应，salmon index常用參數(shù)：
mkdir ~/reference/index/salmon/grcm38
cd ~/reference/index/salmon/grcm38
salmon index -p 12 -t ~/reference/ensembl/grcm38.cdna.fa.gz -i ./
2. salmon定量基因
使用salmon quant命令對clean fastq文件直接進行基因定量，主要參數(shù)如下：
vim 3333_salmon.sh
###################################################
echo -e '\n \n \n ### salmon quant is Working !!! \n \n'
###set##set###set#############
index="/home/reference/index/salmon/grcm38/"
mkdir ~/test/salmon
cd ~/test/salmon
pwd
cat ~/test/idname | while read id
do
echo " '\n !!!!!! Processing sample $id !!!!! '\n"
########single#############################################
# salmon quant -i $index -l A \
# -r ~/test/clean/$(basename $id)_trimmed.fq.gz \
# -p 12 -o $(basename $id)_quant
#######paired#############################################
salmon quant -i $index -l A \
-1 ~/test/clean/${id}_*1*gz \
-2 ~/test/clean/${id}_*2*gz \
-p 12 --output ${id}_quant
###############################################################
done
multiqc ./
echo -e " \n \n \n !!!!ALL WORK DONE !!!!!!!!!!!!!!!!!!!!! \n"
date
nohup bash 3333_salmon.sh >log_333salmon 2>&1 &
運行結(jié)果存放在salmon文件夾下吼虎，里面的quant.sf即為下游分析所需要的文件
三犬钢、獲取基因ID轉(zhuǎn)化的對應文件
由于本次使用的為gencode或ensembl的gtf與cdna文件，因此最后得到的為ensembl_id (gene_id)和 transcript_id思灰，形式為：ENSMUSG00000000001.1 玷犹，而我們下游常用gene symbol進行展示，因此還需要從gtf注釋文件中獲取ensembl_id 洒疚、transcript_id與gene symbol的對應關(guān)系文件箱舞。
方法如下：
vim gtf_geneid2symbol_gencode.sh
#提取gtf注釋文件中g(shù)ene_id等與gene_name的對應關(guān)系,便于下游id轉(zhuǎn)換
gtf="gencode.vM25.chr_patch_hapl_scaff.annotation.gtf"
### gene_id to gene_name
grep 'gene_id' $gtf | awk -F 'gene_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmp
grep 'gene_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmp
paste gene_id_tmp gene_name_tmp >last_tmp
uniq last_tmp >g2s_vm25_gencode.txt
rm *_tmp
### transcript_id to gene_name
grep 'transcript_id' $gtf | awk -F 'transcript_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmp
grep 'transcript_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmp
paste gene_id_tmp gene_name_tmp >last_tmp
uniq last_tmp >t2s_vm25_gencode.txt
rm *_tmp
bash gtf_geneid2symbol_gencode.sh
所得文件如下所示
上游流程到此就結(jié)束了，將最后得到的counts文件夾與g2s_vm25_gencode.txt 或 salmon文件夾與t2s_vm25_gencode.txt下載到本地就可以愉快地進行下游分析了
參考資料
39個轉(zhuǎn)錄組分析工具拳亿，120種組合評估
(https://www.nature.com/articles/s41467-017-00050-4)
轉(zhuǎn)錄組定量可以用替換featureCounts代替HTSeq-count (qq.com)
本實戰(zhàn)教程基于以下生信技能樹分享的視頻：
【生信技能樹】轉(zhuǎn)錄組測序數(shù)據(jù)分析_嗶哩嗶哩_bilibili
【生信技能樹】GEO數(shù)據(jù)庫挖掘_嗶哩嗶哩_bilibili
RNA-seq實戰(zhàn)系列文章：
RNA-seq入門實戰(zhàn)（零）：RNA-seq流程前的準備——Linux與R的環(huán)境創(chuàng)建
RNA-seq入門實戰(zhàn)（一）：上游數(shù)據(jù)下載、格式轉(zhuǎn)化和質(zhì)控清洗
RNA-seq入門實戰(zhàn)（二）：上游數(shù)據(jù)的比對計數(shù)——Hisat2+ featureCounts 與 Salmon
RNA-seq入門實戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣
RNA-seq入門實戰(zhàn)（四）：差異分析前的準備——數(shù)據(jù)檢查
RNA-seq入門實戰(zhàn)（五）：差異分析——DESeq2 edgeR limma的使用與比較
RNA-seq入門實戰(zhàn)（六）：GO愿伴、KEGG富集分析與enrichplot超全可視化攻略
RNA-seq入門實戰(zhàn)（七）：GSEA——基因集富集分析
RNA-seq入門實戰(zhàn)（八）：GSVA——基因集變異分析
RNA-seq入門實戰(zhàn)（九）：PPI蛋白互作網(wǎng)絡構(gòu)建（上）——STRING數(shù)據(jù)庫的使用
RNA-seq入門實戰(zhàn)（十）：PPI蛋白互作網(wǎng)絡構(gòu)建（下）——Cytoscape軟件的使用
RNA-seq入門實戰(zhàn)（十一）：WGCNA加權(quán)基因共表達網(wǎng)絡分析——關(guān)聯(lián)基因模塊與表型
