前言
实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等。生信技能树
htseq-count使用介绍
Usage：htseq-count [options] <sam_file> <gff_file>
主要参数：
-m 计数模型选择，默认union；还有intersection-strict、 intersection-nonempty。
-s 是否链特异性建库，默认yes
-r 排序方法，推荐按name排序，提高效率
-i 指定read计数单位，采用Ensembl GTF时，默认按gene_id计数
-f 输入文件类型 -a 指定mapping质量值过滤read，低的被过滤掉
1503932002082.png
实验操作
ls *.Nsort.bam |while read id;
do
# read计数 -s no -r name
nohup samtools view $id |htseq-count -f sam -s no -r name -i gene_name - ../ref/gencode.vM13.annotation.gtf 1>${id%%.*}.geneCounts 2>${id%%.*}.HTseq.log &;
done
# 合并各个样本的counts值组成表达矩阵
awk 'BEGIN {OFS="\t"}
NR==FNR{
if(FNR==1){f="Gene\t"FILENAME};
a[FNR]=$1;gene[$1]=$2;next}
{if(FNR==1){f=f"\t"FILENAME};
gene[$1]=gene[$1]"\t"$2}
END{print f;for (i in a) print a[i]"\t"gene[a[i]]}'
*.geneCounts > Akap95_expr.txt
ps: htseq-count目前只能单线程跑，多个文件通过循环并行。借助GNU parallel实现htseq-count单文件的并行版。(存在的问题：个别基因counts数目比官方的htseq-count多1个count。原因：对bam分割后，部分reads的两个mates分割进了两个文件，htseq-count把这两个mates都记为了1)
# 借助 GNU parallel，实现 htseq-count 并行版
# -j指定线程；--block文件分块大小(按1G大小分割)；可根据具体配置更改
samtools view Akap95_1.Nsort.bam |parallel -j6 --block 1G --pipe -k "htseq-count -f sam -s no -r name -i gene_name - ../ref/gencode.vM13.annotation.gtf >Akap95_1_p{#}.geneCounts" # time ~30min
# 每个job输出一个文件，互不干扰
# 用awk合并最终总的counts数目， !!!个别bam文件分割处的gene，其counts可能有 1 count的差异(相比普通htseq-count)
awk 'BEGIN {OFS="\t"}NR==FNR{a[FNR]=$1;gene[$1]=$2;next}{gene[$1]+=$2} END{for (i in a) print a[i]"\t"gene[a[i]]}' Akap95_1_p*.geneCounts >Akap95_1_p.geneCounts
如下分别是分割后第2个文件最后一行和第3个文件首行：
1503987659492.png
1503987709136.png
这两个reads都比对到了Aggf1的第5号外显子上
1503988186481.png
1503988844150.png
载入IGV看看：(实验组Akap8(Akap95)表达更高了？？？)
1503992852945.png
GAPDH、β-actin表达情况
1503992170600.png
1503992286983.png
