RNA-seq入门实战（二）：上游数据的比对计数——Hisat2

2022-05-10 本文已影响0人嘿嘿嘿嘿哈

本节概览：

hisat2 + featureCounts:
获取hisat2索引文件，hisat2比对和samtools格式转化，featureCounts计数得到counts表达矩阵

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

获取ensembl_id或transcript_id转化的对应文件

承接上节RNA-seq入门实战（一）
本节介绍使用hisat2或salmon这两种方法进行转录组上游数据的比对和计数。
39个转录组分析工具，120种组合评估
(https://www.nature.com/articles/s41467-017-00050-4)
表明基于hisat2或salmon进行转录本定量都比较优秀。

[Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis | Nature Communications]

一、hisat2 + featureCounts

1. 获取hisat2比对索引文件

index官网下载地址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转化格式

先用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/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 &

比对结果如下：

3. featureCounts 计数得到counts表达矩阵

计数首先要获取gtf注释文件，注意要和hisat2的index文件的基因组版本相对应，如本次为mm10，则gtf文件也必须为mm10或grcm38。
研究人和鼠推荐用gencode数据库的文件GENCODE ，比较常用的还有UCSC的refGene.gtf文件，下载地址在https://hgdownload.soe.ucsc.edu/ （若想下载其他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###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
#######################################################################################            
###生成网页版统计情况
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/grcm38
cd ~/reference/index/salmon/grcm38
salmon 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/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 &

运行结果存放在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="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

所得文件如下所示

上游流程到此就结束了，将最后得到的counts文件夹与g2s_vm25_gencode.txt 或 salmon文件夹与t2s_vm25_gencode.txt下载到本地就可以愉快地进行下游分析了

参考资料
39个转录组分析工具，120种组合评估
(https://www.nature.com/articles/s41467-017-00050-4)
转录组定量可以用替换featureCounts代替HTSeq-count (qq.com)
本实战教程基于以下生信技能树分享的视频：
【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili