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

0.003字数 571阅读 2403

作业要求

1.实现这个功能的软件也很多,还是烦请大家先自己搜索几个教程,入门请统一用htseq-count,对每个样本都会输出一个表达量文件。
2.需要用脚本合并所有的样本为表达矩阵。参考:生信编程直播第四题:多个同样的行列式文件合并起来。
3.对这个表达矩阵可以自己简单在excel或者R里面摸索,求平均值,方差。
看看一些生物学意义特殊的基因表现如何,比如GAPDH,β-ACTIN等等。
来源于生信技能树:http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

实验过程

1.数据准备

文献中测序的是PE,因此我们在对双端数据进行处理时,必须要按reads名进行排序。

# 首先将bam文件按reads名称进行排序(前期是按照默认的染色体位置进行排序的,所以要重新进行排序),这里我们主要以小鼠的数据为例子,不进行人类的测序数据。
$ cd ~/disk2/data/rna-seq/aligned
$ for ((i=59;i<=62;i++));do samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam;done 
# 将排序后的bam文件再次转换成sam格式的文件
$ for ((i=59;i<=62;i++));do samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam;done 

2. reads计数

数据准备已经完成,接下来要使用htseq-count对数据进行reads 计数。
Usage:htseq-count [options] <sam_file> <gff_file>
这里最好使用ensembl的基因组注释文件,小鼠的文件还是需要再次的下载。

#小鼠基因组注释文件的下载,我下载的是m10版本的,与基因组匹配
$ mkdir -p ~/disk2/data/reference/genome/mm10
$ cd ~/disk2/data/reference/genome/mm10
$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz
$ gunzip gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz && rm -rf gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz
# 之前的环境已经全部搭建完成,直接就可以使用htseq-count
$ mkdir -p ~/disk2/data/rna-seq/matrix && cd ~/disk2/data/rna-seq/matrix
$ for ((i=59;i<=62;i++));do htseq-count ~/disk2/data/rna-seq/aligned/SRR35899${i}_count.sam ~/disk2/data/reference/genome/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf;done 

最后得到4个矩阵文件
reads 计数文件

count文件的格式:基因名一列+reads计数一列
count文件结构

这一部分主要参考的是老司机Jimmy大神的博客:http://www.bio-info-trainee.com/244.html

3.合并表达矩阵

使用R将这四个文件进行合并,得到最后的融合表达矩阵,R语言代码:

>options(stringsAsFactors = FALSE) 
# 首先将四个文件分别赋值:control1,control2,rep1,rep2
>control1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
>control2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2")) 
>rep1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951")) 
>rep2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
# 将四个矩阵按照gene_id进行合并,并赋值给raw_count
>raw_count <- merge(merge(control1, contol2, by="gene_id"), merge(rep1,rep2, by="gene_id"))
# 需要将合并的raw_count进行过滤处理,里面有5行需要删除的行,在我们的小鼠的表达矩阵里面,是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter
>raw_count_filt <- raw_count[-48823:-48825, ]
>raw_count_filter <- raw_count_filt[-1:-2, ]
# 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因,只能是ENSMUSG00000024045的整数,没有小数点,所以需要进一步替换为整数的形式。
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空,并赋值给ENSEMBL
>ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt1$gene_id) 
# 将ENSEMBL重新添加到raw_count_filt1矩阵
>row.names(raw_count_filter) <- ENSEMBL
# 看一些基因的表达情况,在UniProt数据库找到AKAP95的id,并从矩阵中找到访问,并赋值给AKAP95变量
>AKAP95 <- raw_count_filter[rownames(raw_count_filt1)=="ENSMUSG00000024045",]
# 查看AKAP95
>AKAP95

  • raw_count_filter结构


    raw_count_filter数据结构
  • AKAP95基因的表达情况,可以看到表达量是提高
AKAP95reads计数情况

这一部分主要参考徐州更同学的R语言代码:http://www.jianshu.com/p/e9742bbf83b9
我的数据是用的小鼠的测序结果,所以我修改了部分的内容,更好的进行。

PS:前段时间一直在忙于实验,没有及时去完成这一部分的内容,今天正好是星期天,所以就打算抽出一点时间来完成这一部分的学习笔记。还有一个原因就是因为自己好像不会了,没法再进行下去了,说到第还是有些害怕了,实在惭愧哈,继续努力加油。

推荐阅读更多精彩内容