转录组学习笔记06-reads计数

在这一部分中，对前面生成的SAM文件进行定量，用到的工具是htseq，如下所示：

mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
    htseq-count -s no -r pos -f bam bio/bio/RAW_DATA/RNA-seq/01/outmice/SRR35899${i}_sorted.bam reference/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > bio/bio/RAW_DATA\RNA-seq/01/reads/SRR35899${i}count 2> RNA-Seq/matrix/SRR35899${i}.log
done

生成的数据如下所示：

mark

用Notepad++打开其中的一个文件，例如SRR3589959count，如下所示：

mark

数据分为两列，第1列是基因名，第2列是reads数。现在就要把这4个文件给合并起来，构成一个矩阵，如下所示

options(stringsAsFactors = FALSE) 
# 首先将四个文件分别赋值：control1，control2，rep1，rep2
control1 <- read.table("F:/RAW_DATA/RNA-seq/01/reads/SRR3589959count", sep="\t", col.names = c("gene_id","control1"))
control2 <- read.table("F:/RAW_DATA/RNA-seq/01/reads/SRR3589961count", sep="\t", col.names = c("gene_id","control2")) 
rep1 <- read.table("F:/RAW_DATA/RNA-seq/01/reads/SRR3589960count", sep="\t", col.names = c("gene_id","akap951")) 
rep2 <- read.table("F:/RAW_DATA/RNA-seq/01/reads/SRR3589962count", sep="\t",col.names = c("gene_id","akap952"))
# 将四个矩阵按照gene_id进行合并，并赋值给raw_count
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(rep1,rep2, by="gene_id"))
head(raw_count)
raw_count_filt <- raw_count[-c(1:5),]# 删除前五行
# 因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id) 
# 将ENSEMBL重新添加到raw_count_filt1矩阵
row.names(raw_count_filt) <- ENSEMBL
# 看一些基因的表达情况，在UniProt数据库找到AKAP95的id，并从矩阵中找到访问，并赋值给AKAP95变量
AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSMUSG00000024045",]
# 查看AKAP95
AKAP95

参考资料