表达定量

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 &