转录组入门（6）：reads计数

2017-08-13 本文已影响737人 lxmic

作业要求

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