转录组入门(6): reads计数

96
宇天_ngs
2017.09.14 11:54* 字数 455

前言

实现这个功能的软件也很多,还是烦请大家先自己搜索几个教程,入门请统一用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
RNA-seq