生信技能树-day18 转录组上游分析-比对、定量

生信菜鸟团

发布于 2024-06-25 20:20:38

940

发布于 2024-06-25 20:20:38

文章被收录于专栏：生信菜鸟团生信菜鸟团

从我们生信技能树历年的几千个马拉松授课学员里面募集了一些优秀的创作者，某种意义来说是传承了我们生信技能树的知识整理和分享的思想！

今天的是三周合计15天的数据挖掘授课学员一点一滴整理的授课知识点笔记哦，还有互动练习题哈，欢迎大家点击文末的阅读原文去关注我们学员的公众号哦！

数据比对

目标：使用两个软件对fq数据进行比对，得到比对文件sam/bam，并探索比对结果。

需要准备：

参考基因组文件fasta
参考基因组注释文件gff/gtf

参考基因组

参考基因组准备:注意参考基因组版本信息，可以用ncbi或者Ensembl数据库，一般用Ensembl数据库，更新较快，基因没有重名

（服务器中已经下载好参考基因组，此处只要了解一下怎么下载即可）

ncbi：https://www.ncbi.nlm.nih.gov/projects/genome/guide/human/index.shtml
Ensembl：http://asia.ensembl.org/index.html

# 具体操作：进入官网，右键复制下载连接，黏贴然后运行对应的脚本
# http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

# 进入到参考基因组目录
mkdir -p $HOME/database/GRCh38.111
cd $HOME/database/GRCh38.111

# 下载基因组序列axel  curl  
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz >dna.log &

# 下载转录组序列
nohup axel -n 100  http://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz >rna.log &

# 下载基因组注释文件
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gtf/homo_sapiens/Homo_sapiens.GRCh38.111.gtf.gz >gtf.log &

nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gff3/homo_sapiens/Homo_sapiens.GRCh38.111.gff3.gz >gff.log&

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
nohup gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz Homo_sapiens.GRCh38.cdna.all.fa.gz >unzip.log &

fasta介绍

• 以“>”开头，序列名称&序列描述 • 序列中允许空格，换行，空行，直到下一个“>”，表示该序列结束

gff/gtf文件介绍

Generic Feature Format，主要用来描述基因的结构与功能信息，对基因组进行注释，目前多用版本为GFF3

格式：文本文件，共9列

第九列的详解

GTF文件

gene transfer format，主要是用来对基因进行注释，前八个字段与GFF相同（有一些小的差别），重点在第九列的不同

Ensembl数据库的ID号

格式：

ENS【species prefix】【feature type prefix】【a unique eleven digit number】

mouse gene：ENSMUSG###########

课后作业

1.fastq与fasta文件转换

# 答案1
zless -S SRR1039511_1_val_1.fq.gz |awk '{ if(NR%4==1){print">" substr($0,2)} if(NR%4==2){print}  }' | less -S

# 答案2
zless -S SRR1039510_1_val_1.fq.gz |paste - - - - |cut -f 1,2 |tr '@' '>' |tr '\t' '\n' |less -S

2.从gff或者gft文件中获取基因的ID与symbol对应关系，以及biotype类型

应用：ID与symbol转换本地化，不依赖于第三方工具和软件包，并可以根据biotype类型区分mRNA，lncRNA以及miRNA等信息。

# 从gff或者gft文件中获取ID与symbol对应关系，以及biotype类型
zless -S Homo_sapiens.GRCh38.104.chr.gtf.gz  |awk -F'\t' '{if($3=="gene"){print$9}}' |awk -F';' '{print$1,$3,$5}' |awk '{print$2"\t"$4"\t"$6}' |sed 's/"//g' |grep 'protein_coding' >protein_coding_id2name.xls

数据比对的过程

建索引：为了将短片段快速比对到基因组上的某一个位置
比对参考基因组，结果生成sam文件
sam转bam
bam建索引

比对:hisat2

hisat2的主要参数

其中链特异性参数和所测的rna是什么类型有关，可以咨询公司

## ----1.构建索引
# 进入参考基因组目录
cd $HOME/database/GRCh38.105
# Hisat2构建索引
# vim Hisat2Index.sh
mkdir Hisat2Index
fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
hisat2-build -p 5 -f ${fa} Hisat2Index/${fa_baseName}

# 运行
nohup bash Hisat2Index.sh >Hisat2Index.sh.log &

## 构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
## 后续索引可直接使用服务器上已经构建好的进行练习

## ----2.比对
# 进入比对文件夹
cd $HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2

## 单个样本比对，步骤分解
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/
outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2

hisat2 -p 5  -x  ${index} \
    -1 ${inputdir}/SRR1039510_1_val_1.fq.gz \
       -2 ${inputdir}/SRR1039510_2_val_2.fq.gz \
       -S ${outdir}/SRR1039510.Hisat_aln.sam
# 其中-p 5是指5个线程，-x index是指比对的参考基因组的路径，-1和-2是指输入的cleandata的read1和read2，-S outdir是指生成的sam文件

# 98.42% overall alignment rate 指总比对率，这个指标非常重要，表示有多少可以比对上参考基因组，一般需要大于80%，但根据样本类型不同，比如外泌体一般只有40-50%左右

# 22844 (92.94%) aligned concordantly exactly 1 time 指唯一比对率，越高越好

## ----3.sam转bam
samtools sort --threads 5 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

## ----4.对bam建索引
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai


# 多个样本批量进行比对，排序，建索引
# Hisat.sh内容： 注意命令中的-，表示占位符，表示|管道符前面的输出。
## 此处索引直接使用服务器上已经构建好的进行练习
# vim Hisat.sh
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/
outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2

cat ../../data/cleandata/trim_galore/ID | while read id
do
hisat2 -p 5 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log  | samtools sort -@ 3 -o ${outdir}/${id}.Hisat_aln.sorted.bam - &&  samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai
done

# ${outdir}/${id}.Hisat_aln.sorted.bam - &&  samtools index
# 以上命令中的-指占位符，表示前一个任务的输出结果通过管道符传递给后一个命令，并指定位置，&&指多个命令串联，只有前一个命令运行成功后才会运行后面的命令

# 提交后台运行
nohup bash Hisat.sh >Hisat.log &

# 统计比对情况
multiqc -o ./ SRR*log

比对结果文件SAM的解释

SAM(The Sequence Alignment/Map format)格式，即序列比对文件格式，详细介绍见：http://samtools.github.io/hts-specs/SAMv1.pdf

FLAG详解：http://broadinstitute.github.io/picard/explain-flags.html FLAG: Combination of bitwise FLAGs，表示比对的详细情况

想知道某一个值是什么意思的话可以直接去网页上搜索对应的值

CIGAR详解：CIGAR string，简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果

扩展内容：sam/bam应用

samtools工具：http://www.htslib.org/doc/samtools.html
Samtools常用命令的总结：https://www.bioinfo-scrounger.com/archives/245/ https://www.cnblogs.com/xiaofeiIDO/p/6805373.html

samtools工具使用

##----view查看bam文件
samtools view SRR1039510.Hisat_aln.sorted.bam
samtools view -H SRR1039510.Hisat_aln.sorted.bam
samtools view -h SRR1039510.Hisat_aln.sorted.bam

##----index对bam文件建索引
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai

##----flagstat统计比对结果
samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam

##----sort排序 sam转bam并排序
samtools sort -@ 3 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

##----depth统计测序深度
# 得到的结果中，一共有3列以指标分隔符分隔的数据，第一列为染色体名称，第二列为位点，第三列为覆盖深度
samtools depth SRR1039510.Hisat_aln.sorted.bam >SRR1039510.Hisat_aln.sorted.bam.depth.txt

##----计算某一个基因的测序深度
# 找到基因的坐标
zless -S Homo_sapiens.GRCh38.95.gff3.gz |awk '{if($3=="gene")print}' |grep 'ID=gene:ENSG00000186092' |awk '{print $1"\t"$4"\t"$5}' >ENSG00000186092.bed
samtools depth  -b  ENSG00000186092.bed SRR1039510.Hisat_aln.sorted.bam >ENSG00000186092.bed.depth

# 如何找到多比对的reads，flag值的理解
# (0x100) 代表着多比对情况，所以直接用samtools view -f 0x100可以提取 multiple比对的 情况

比对:subjunc

参考文档：http://subread.sourceforge.net/

常用参数

## ----1.构建索引
# 进入参考基因组目录
cd $HOME/database/GRCh38.105
# subjunc构建索引，构建索引时间比较长大约40分钟左右，建议写成sh脚本提交后台运行
## 后续索引可直接使用服务器上已经构建好的进行练习
# vim SubjuncIndex.sh
mkdir SubjuncIndex
fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
subread-buildindex -o SubjuncIndex/${fa_baseName} ${fa}

# 运行
nohup bash SubjuncIndex.sh >SubjuncIndex.sh.log &


## ----比对
# 进入文件夹
cd $HOME/project/Human-16-Asthma-Trans/Mapping/Subjunc

# vim subjunc.sh
index=/home/t_rna/database/GRCh38.104/SubjuncIndex/GRCh38.dna
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Subjunc

cat ../../data/cleandata/trim_galore/ID | while read id
do
  subjunc -T 10 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R ${inputdir}/${id}_2_val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 -o ${outdir}/${id}.Subjunc.sorted.bam  ${outdir}/${id}.Subjunc.bam && samtools index ${outdir}/${id}.Subjunc.sorted.bam ${outdir}/${id}.Subjunc.sorted.bam.bai
done

# 运行
nohup bash subjunc.sh >subjunc.log &
# 运行结果中每个样本（read1+2）会生成一个subjunc.log
# .indel.vcf是插入（in）和缺失（del）的信息

表达定量

featureCounts

使用featureCounts对bam文件进行定量。

官网：http://bioinf.wehi.edu.au/featureCounts/

featureCounts 需要两个输入文件：

比对产生的BAM/ SAM文件
区间注释文件（GTF格式，SAF格式）

cd $HOME/project/Human-16-Asthma-Trans/Expression/featureCounts

## 定义输入输出文件夹
gtf=/home/t_rna/database/GRCh38.104/Homo_sapiens.GRCh38.104.chr.gtf.gz
inputdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2/

# featureCounts对bam文件进行计数
featureCounts -T 6 -p --countReadPairs -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.sorted.bam
# -T 6为6个线程，--countReadPairs -t exon是指想要统计的是exon，-g gene_id为统计的meta-feature为gene_id，-a $gtf指输入的参考基因组注释文件，-o all.id.txt指输出文件，最后跟输入文件

# 对定量结果质控
multiqc all.id.txt.summary

featureCounts的结果解析

# 得到表达矩阵txt文件，需要进一步处理为行为基因，列为样本的原始表达矩阵raw counts
# 处理表头，/home/t_rna/要换成自己的路径
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#/home/t_rna/project/Human-16-Asthma-Trans/Mapping/Hisat2//##g' |sed 's#.Hisat_aln.sorted.bam##g' >raw_counts.txt
# less出表达矩阵，然后grep将表头注释行去掉，cut取出第一列和第7列及以后的列，sed用连续的三个相同字符（因为/太多了此处不用/）使用命令s/pattern/new/[flags]替换字符串，即将/home/t_rna/project/Human-16-Asthma-Trans/Mapping/Hisat2//替换为空，g表示处理每一行，然后将结果又传递给sed，将.Hisat_aln.sorted.bam替换为空，最后将结果写入raw_counts.txt

# sed可以用任意连续三个相同字符分隔，比如：
sed s///
sed s###
sed s%%%

# 列对齐显示
head raw_counts.txt  |column -t

salmon定量

Salmon可以快速从fastq快速得到基因表达，号称不用比对，直接定量

Salmon参考文档：https://salmon.readthedocs.io/en/latest/

-t：参考基因组fasta文件，可以接受压缩格式

-i：存储索引的文件夹名

##----构建索引
## 后续索引可直接使用服务器上已经构建好的进行练习
cd $HOME/database/GRCh38.105
nohup salmon index -t Homo_sapiens.GRCh38.cdna.all.fa -i salmon_index >salmon-index.log &

cd $HOME/project/Human-16-Asthma-Trans/Expression/Salmon

##----运行
# 编写脚本，使用salmon批量对目录下所有fastq文件进行定量
# vim salmon.sh
index=/home/t_rna/database/GRCh38.104/salmon_index
input=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
outdir=$HOME/project/Human-16-Asthma-Trans/Expression/Salmon
cat ../../data/cleandata/trim_galore/ID |while read id 
do
  salmon quant -i ${index} -l A -1 ${input}/${id}_1_val_1.fq.gz -2 ${input}/${id}_2_val_2.fq.gz -p 5 -o ${outdir}/${id}.quant
done

salmon quantmerge --quants {SRR1039510.quant,SRR1039511.quant,SRR1039512.quant} --names {SRR1039510,SRR1039511,SRR1039512} --column numreads  --output raw_count.txt

# tpm矩阵
salmon quantmerge --quants {SRR1039510.quant,SRR1039511.quant,SRR1039512.quant} --names {SRR1039510,SRR1039511,SRR1039512} --column tpm  --output tpm.txt

# 以上所有命令可以直接放进salmon.log中，然后直接：
# 后台运行脚本
nohup bash salmon.sh >salmon.log &

# 其中SRR1039510.quant,SRR1039511.quant,SRR1039512.quant如果少量还可以自己写，如果多的话怎么生成？
##----合并表达矩阵
# 原始count值矩阵
# --quants：ls -d *quant |tr '\n' ',' |sed 's/,$//' |awk '{print "{" $0 "}" }'
# --quants：ls -d *quant |sed ':a;N;s/\n/,/;t a;'|awk '{print "{"$0"}"}'
# --quants：ls -d *quant |xargs  |tr ' ' ',' |awk '{print "{"$0"}"}'
# --names：ls -d *quant |tr '\n' ',' |sed 's/,$//' |awk '{print "{"$0"}"}' |sed 's/.quant//g'

## 生成完整版命令：ls -d *quant |tr '\n' ',' |sed 's/,$//' |awk '{print "{"$0"}"}'  | perl -ne 'chomp;$a=$_;$a=~s/\.quant//g;print"salmon quantmerge --quants $_ --names $a --column numreads  --output raw_count.txt \n";'