RNA-seq(4) :采用Feature counts进行reads计数，合并矩阵,并进行基因ID注释-学习笔记

这段时间去读了一些人文类的书籍，没更新，一看到快到八月份，继续学点分析方面的东西。

主要参考网站：
RNA-seq(6): reads计数，合并矩阵并进行注释
原创10000+生信教程大神给你的RNA实战视频演练
Biomat对基因ID转换
R语言删除/添加数据框中的某一行/列
R-按照行名或列名删除对应的行列

序列对比之后，按照以往的分析程序，接下来要看reads数目多少了。最常使用的软件是htseq。可以参考用htseq-count对reads计数并合并矩阵。但是这个方法真的很费时间，所以找到了一个便捷的工具，Featurecounts。

还是要先下载，因为之前我在学习RNA-seq的时候下载了htseq,但是没有下载这个。

#featureCounts 在 subread 这个包里面
conda install subread
#建立一个软链接,方便直接调用
 ln -s ~/miniconda3/bin/featureCounts ~/.local/bin 

需要说明的是，一定要把bam文件按照顺序（默认是name）排序，代码在上一个学习笔记里有，但是为了方便，继续抄下来。

 for ((i=22;i<=23;i++));do samtools sort SRR53377${i}.bam -o SRR53377${i}
.sorted_bam;done
#随后用featurecounts计数
for ((i=19;i<=23;i++));do featureCounts -T 5 -p -t exon -g gene_id  -a /mnt/d/biotech/chip/ref/gencode.vM10.annotation.gtf.gz -o ${i}.counts.txt SRR53377${i}_sorted.bam;done

下面是运行的结果：
4min就比对了一个counts文件，要是htseq真的不好说，参考链接上说的是5h,4个文件这么久，真的害怕极了。还好有这个小神器，哈哈。

5个counts 文件大概用了40min。

接下来，看一下文件可以去文件夹里用记事本打开，也可以用代码。

 tail -n 5 19.counts.txt
 head -n 6 22.counts.txt

结果分别如下：


这儿就会神奇的发现，没有htseq计数之后产生的那些需要删除的东西，省了很多的事情，比如这种：

接下来，还可以查看一下总的结果

wc -l *.counts.txt

结果如下

下面的操作都要在Rstudio里进行了。

  #更改一下当前的文件路径设置
  setwd("D:/biotech/RNA/aligned") 
getwd()   #查看文件路径
#加载数据，19——21是实验组，22-23是对照组
options(stringsAsFactors = FALSE)
control1<-read.table("22.counts.txt",sep = "\t",header =T)
#这个操作是把下面的数据集第七列的名称SRR537722改成control1
colnames(control1)[7]<-("contorl1")
#按照上边两个代码加载其他的数据集这儿省略

#删除一些不太重要的信息比如基因所在的位置和长度等信息，即数据集第2-6列，第一列是GNEE ID,第7列是样品名，比如control2，或者treat3
treat1 <- treat1[,-(2:6)]

#合并上述五个数据集，但是merge这个函数只可以合并两个数据集，所以采用了下述代码
raw_count<-merge(control1,control2,by="Geneid")
raw_count2<-merge(treat1,treat2,by="Geneid")
raw_count2<-merge(raw_count2,treat3,by="Geneid")
raw_count <- merge(raw_count, raw_count2, by="Geneid")
#也可以按照参考网址上的代码修改
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))

#merge之后顺序会变，删除的时候看下删除的是哪些行，这个是htseq后续的操作，featurecounts用不到，在这儿只是记一下R的数据集如何删除行或列
raw_count_filt <- raw_count[-1:-5,]

RNA-seq(6): reads计数，合并矩阵并进行注释 为 #参考链接

#由于数据库里的基因ID没有小数点，所以需要进一步替换为整数的形式。
#将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count$Geneid) 
#参考链接里是把 ENSEMBL 作为列名加入数据集，这儿做了改动，使ENSEMBL作为一列的内容插入到数据集第2列中，以便后续需要时进行更改，再把数据集第一列，也就是带有小数点之后数字的一列去掉。
raw_count <- cbind(raw_count[,1],ENSEMBL,raw_count[,2:ncol(raw_count)])
raw_count<-raw_count[,-1]

（可选） 对基因进行注释，其实这一步可以在这个地方省略，这一步在这儿做的目的是提前对比一下对照组和实验组的reads读数差别，提前进行一个验证，做差异分析之后还需要进行一次注释）

library('biomaRt')
#数据集的第一列的ID名为参照去搜索对应的基因名称
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
my_ensembl_gene_id<-raw_count[,1]
mms_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name'),filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
#为了把基因的名称和数据集进行合并，要把两个数据集其中一列的名称改成一致的，比如都改成  ensembl_gene_id
colnames(raw_count)[1]<-("ensembl_gene_id")
readout=merge(mms_symbols,raw_count,by="ensembl_gene_id")
#查看一个基因的表达情况
Akap8 <- readout[readout$external_gene_name=="Akap8",]
Akap8
#结果如下，很明显处理组表达增高
ensembl_gene_id external_gene_name contorl1 contorl2 treat1 treat2 treat3
4261 ENSMUSG00000024045              Akap8     1194     1002   1636   1639   1496

#如果需要的话，记得保存
write.csv(raw_count,"readout.csv",row.names=FALSE)