生信技能樹知識整理實習生招募，讓我走大運結識了幾位優秀小夥伴！大家開始根據我的ngs組學視頻進行一系列公共數據集分析實戰，其中幾個小夥伴讓我非常驚喜，不需要怎麼溝通和指導，就默默的完成了一個實戰！

他前面的分享是：

• Counts FPKM RPKM TPM CPM 的轉化

• 獲取基因有效長度的N種方

下面是他對我們b站轉錄組視頻課程的詳細筆記

本節概覽：

1. hisat2 + featureCounts:
   獲取hisat2索引文件，hisat2比對和samtools格式轉化，featureCounts計數得到counts表達矩陣

2. Salmon:
   salmon index 用cdna.fa建立索引，salmon quant對clean fastq文件直接進行基因定量

3. 獲取ensembl_id或transcript_id轉化的對應文件

承接上節RNA-seq入門實戰（一）本節介紹使用hisat2或salmon這兩種方法進行轉錄組上遊數據的比對和計數。39個轉錄組分析工具，120種組合評估(/articles/s41467-017-00050-4)表明基于hisat2或salmon進行轉錄本定量都比較優秀。

一、hisat2 + featureCounts

1. 獲取hisat2比對索引文件

index官網下載地址Download | HISAT2 (daehwankimlab.github.io)，下載并解壓所需的 mm10 或 grcm38 的index文件

mkdir -p ~/reference/index/hisat/cd ~/reference/index/hisat/wget /hisat/mm10_genome.tar.gztar -zxvf *tar.gz rm *tar.gz

2. hisat2比對和samtools轉化格式

先用hisat2比對基因組得到sam文件，再用samtools sort将sam文件格式轉化與排序為bam文件（bam相當于二進制版的sam），之後samtools index建立索引（用于後續IGV内可視化），最後samtools flag 統計文件比對情況保存在文本中。其中samtools index與samtools flag為非必須步驟，可略過。sam相當于是中間文件比較占存儲空間，可在轉化為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/flagcd ~/test/align/pwdcat ~/test/idname | while read iddo 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}.samdone
#### samtools index and flagstat ####echo -e " \n \n \n samtools index and flagstat \n " cd -p ~/test/align/flagpwd#ls ~/test/align/*.bam | xargs -i samtools index -@ 12 {} ls ~/test/align/*.bam | while read id ;\do samtools flagstat -@ 12 $id > $(basename $id '.bam').flagstatdonemultiqc ./
echo -e " \n \n \n 333# ALL Work Done !!!\n \n \n "date

nohup bash 3_align2sam2bam_hisat2.sh >log_3 2>&1 &

比對結果如下：

3. featureCounts 計數得到counts表達矩陣

計數首先要獲取gtf注釋文件，注意要和hisat2的index文件的基因組版本相對應，如本次為mm10，則gtf文件也必須為mm10或grcm38。
研究人和鼠推薦用gencode數據庫的文件GENCODE，比較常用的還有UCSC的refGene.gtf文件，下載地址在/（若想下載其他gtf文件則将網址中mm10替換即可，如hg38）。
featureCounts詳細使用方法見轉錄組定量可以用替換featureCounts代替HTSeq-count (qq.com)，常用參數如下：

代碼如下：

vim 4_counts_featurecounts.sh

###########################################echo -e "\n \n \n444# Count #featureCounts !!! \n \n \n"date#####set####set###setgtf='/home/reference/gtf/gencode/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf'#gtf='/home/reference/gtf/UCSC/mm10.refGene.gtf.gz'
mkdir -p ~/test/countscd ~/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####################################################################################### ###生成網頁版統計情況multiqc ./
echo -e " \n \n \n ALL WORK DONE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \n "date

nohup bash 4_counts_featurecounts.sh >log_4 2>&1 &

查看計數的統計情況，匹配率在71-80%左右，還可以

運行結果保存在counts文件夾下，裡面的counts.txt就是我們下遊分析所需要的文件啦

---

二、Salmon——直接對基因進行定量的工具

與hisat2不同，salmon不經過比對計數步驟而是直接對基因進行定量，如果不研究新轉錄本，用salmon方法可以更快更方便得到基因定量信息。

1. 建立salmon索引

先下載參考轉錄本序列cDNA.fa文件，在ensembl官網選擇相應文件 Index of /pub/release-102/fasta/mus_musculus/cdna/ (ensembl.org)，使用salmon index 建立索引文件，salmon index常用參數：

mkdir ~/reference/index/salmon/grcm38cd ~/reference/index/salmon/grcm38salmon index -p 12 -t ~/reference/ensembl/grcm38.cdna.fa.gz -i ./

2. salmon定量基因

使用salmon quant命令對clean fastq文件直接進行基因定量，主要參數如下：

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/salmoncd ~/test/salmonpwdcat ~/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 &

運行結果存放在salmon文件夾下，裡面的quant.sf即為下遊分析所需要的文件

---

三、獲取基因ID轉化的對應文件

由于本次使用的為gencode或ensembl的gtf與cdna文件，因此最後得到的為ensembl_id (gene_id)和 transcript_id，形式為：ENSMUSG00000000001.1 ，而我們下遊常用gene symbol進行展示，因此還需要從gtf注釋文件中獲取ensembl_id 、transcript_id與gene symbol的對應關系文件。
方法如下：

vim gtf_geneid2symbol_gencode.sh

#提取gtf注釋文件中gene_id等與gene_name的對應關系,便于下遊id轉換

#提取gtf注釋文件中gene_id等與gene_name的對應關系,便于下遊id轉換gtf="gencode.vM25.chr_patch_hapl_scaff.annotation.gtf"
### gene_id to gene_namegrep 'gene_id' $gtf | awk -F 'gene_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmpgrep 'gene_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmppaste gene_id_tmp gene_name_tmp >last_tmpuniq last_tmp >g2s_vm25_gencode.txtrm *_tmp
### transcript_id to gene_namegrep 'transcript_id' $gtf | awk -F 'transcript_id \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_id_tmpgrep 'transcript_id' $gtf | awk -F 'gene_name \"' '{print $2}' |awk -F '\"' '{print $1}' >gene_name_tmppaste gene_id_tmp gene_name_tmp >last_tmpuniq last_tmp >t2s_vm25_gencode.txtrm *_tmp

bash gtf_geneid2symbol_gencode.sh

所得文件如下所示:

---

上遊流程到此就結束了，将最後得到的counts文件夾與g2s_vm25_gencode.txt 或 salmon文件夾與t2s_vm25_gencode.txt下載到本地就可以愉快地進行下遊分析了

參考資料

39個轉錄組分析工具，120種組合評估(/articles/s41467-017-00050-4)

轉錄組定量可以用替換featureCounts代替HTSeq-count (qq.com)

本實戰教程基于以下生信技能樹分享的視頻：

【生信技能樹】轉錄組測序數據分析_哔哩哔哩_bilibili

【生信技能樹】GEO數據庫挖掘_哔哩哔哩_bilibili

文末友情宣傳

強烈建議你推薦給身邊的博士後以及年輕生物學PI，多一點數據認知，讓他們的科研上一個台階：

• 數據挖掘（GEO,TCGA,單細胞）2022年5~6月場，快速了解一些生物信息學應用圖表
• 生信入門課-2022年5~6月場，你的生物信息學第一課