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文件进行计数
# 如果是Version 2.0.1,使用下面的代码
featureCounts -T 6 -p -t exon -g gene_id -a $gtf -o all.count.txt $inputdir/*.sorted.bam
#如果是Version v2.0.3,使用下面的代码:
featureCounts -T 6 -p --countReadPairs -t exon -g gene_id -a $gtf -o all.count.txt $inputdir/*.sorted.bam
#-T线程;-p --countReadPairs针对双端测序数据;-g 在基因水平进行定量
-a 注释文件 -o 输出的定量后的文件 后面是用到的bam文件
# 生成结果 图1;结果解析图2
#对定量结果质控
multiqc all.count.txt.summary
# 得到表达矩阵
从图1到图3
# 处理表头,/home/t_rna/要换成自己的路径
grep -v '#' all.count.txt |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
sed s///
sed s###
sed s%%%
# 列对齐显示
head raw_counts.txt |column -t
##----构建索引
## 后续索引可直接使用服务器上已经构建好的进行练习
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
如图4
##----合并表达矩阵
# 原始count值矩阵 (替换脚本中的SRR....*.quant 这些文件的名称)
# --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";'
#按照column numreads 合并
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
# 后台运行脚本
nohup sh salmon.sh >salmon.log &
需要的结果文件是quant.sh 图5(之后需要合并)
----来自生信技能树----
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。