生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的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月场，你的生物信息学第一课