对测序得到的reads进行计数，即基因表达的定量过程。根据reads和基因位置的overlap，以此来判断reads到底属于哪一个基因，同时对该reads总数进行计数，生成counts矩阵。
今日内容
1.HTSeq-count对reads进行计数
2.R语言完成counts矩阵的合并
1. HTSeq–count对reads进行计数
首先了解HTseq用法，参数说明如下：
usage: htseq-count [options] alignment_file gff_file
positional arguments:samfilenames Path to the SAM/BAM files containing the mapped reads.If '-' is selected, read from standard inputfeaturesfilename Path to the file containing the featuresoptional arguments:-h, --help show this help message and exit-f {sam,bam}, --format {sam,bam}type of <alignment_file> data, either 'sam' or 'bam'(default: sam) 输入文件类型sam/bam；默认为sam文件 -r {pos,name}, --order {pos,name}'pos' or 'name'. Sorting order of <alignment_file>(default: name). Paired-end sequencing data must besorted either by position or by read name, and thesorting order must be specified. Ignored for single-end data.输入文件的排序方式，默认按read名排序--max-reads-in-buffer MAX_BUFFER_SIZEWhen <alignment_file> is paired end sorted byposition, allow only so many reads to stay in memoryuntil the mates are found (raising this number willuse more memory). Has no effect for single end orpaired end sorted by name-s {yes,no,reverse}, --stranded {yes,no,reverse}whether the data is from a strand-specific assay.Specify 'yes', 'no', or 'reverse' (default: yes).'reverse' means 'yes' with reversed strandinterpretation-a MINAQUAL, --minaqual MINAQUALskip all reads with alignment quality lower than thegiven minimum value (default 10）剔除mapping quality值低于阈值的read-t FEATURETYPE, --type FEATURETYPEfeature type (3rd column in GFF file) to be used, allfeatures of other type are ignored (default, suitablefor Ensembl GTF files: exon)前一期利用SAMtools工具对比对产生的sam文件进行格式转换并排序，得到了按染色体位置进行排序的bam文件${id}.sorted.bam。
因此输入文件格式-f bam，文件名${id}.sorted.bam ;排序方式-r pos ;文献中提到reads计数时剔除了mapping quality < 30的低质量reads，即-a 30。
输入bam文件路径：./aligned/${id}.sorted.bam
参考基因组gtf索引文件./genome_gtf/gencode.v41lift3_.gtf
生成的count文件存放于./aligned/
for id in `seq 56 62`;
do
#bsub -J ${id} -n 8 -o %J.${id}.out -e %J.${id}.err -R span[hosts=1] -q normal \
htseq-count -r pos -a 30 \
-f bam ./aligned/${id}.sorted.bam \
./genome_gtf/gencode.v41lift37.basic.annotation.gtf >./aligned/${id}.count
done运行结束将产生.count结尾的文件，head查看一下count文件前15行，了解一下数据结构。结构是一个二维矩阵，第一列为Ensemble ID，小数点后面的数字部分为ID版本信息；第二列即为reads数，可理解为表达量。
至此基于Linux系统的RNA-seq数据上游分析流程基本完成；下游分析主要是可视化过程，依赖R语言来实现～
2.基于Rstudio的下游分析——counts矩阵的准备
回观count文件的数据结构，你会发现我们面临一个问题：第一列ENS编号并不是我们所悉的基因名称(gene_id)。因此在下游分析前还需完成Ensemble ID版本号的去除及向Gene_ID的转换，如果是含有多个样本的count矩阵，还需将多个count矩阵合并为一个多维矩阵。以上流程可通过Rstudio来完成，使用的R包包括stringr（str_split ）、biomaRt(getBM)、merge函数。
2.1. 软件准备
R包使用逻辑：安装——R包调用——函数调用；对于没有安装的R包可通过'install.package("package name")'命令进行下载安装。但使用以上命令安装biomaRt及其依赖包AnnotationDbi时将报错，无法获取安装包信息，可使用BioMannager进行安装，亲测有效
install.package('stringr‘）
BiocManager::install("“biomaRt”")
BiocManager::install("AnnotationDbi")#软件包下载
library(biomaRt)#R包调用
library(stringr)2.2.多个count矩阵的合并（merge）
由于这个系列所用测序数据包含7个测序数据（SRR56～62），在reads计数时会各自生成一个count矩阵，下游分析时一个一个处理count文件比较麻烦，需将七个count文件进行合并；由于reads计数时使用的是同一个基因组索引文件，因此生成的count文件Ensemble ID必定是一致的，这就方便我们以第一列ID为交集索引实现文件合并。使用的函数为基础函数merge。
上图为SRR3589960～61 count文件部分行名展示
注⚠️：为节省时间，最初只选用了SRR3589960～SRR3589962三个mRNA测序数据进行测试，因此之后的流程只选择这3个测序数据对应的count矩阵进行演示；
首先通过help（'merge')查看函数基本参数用法
#用法Usage:
merge(x, y, by = intersect(names(x), names(y)),by.x = by, by.y = by, all = FALSE, all.x = all, all.y = all,sort = TRUE, suffixes = c(".x",".y"), no.dups = TRUE,incomparables = NULL, ...)
#x，文件1；y，文件2；by参数对于文件1与文件2的交集列（行）名
#all，如果为 TRUE，x中的每一行在y中没有匹配的行将被添加到输出结果中且这些行通常用NA填充；默认为 FALSE#根据count数据结构，确定使用参数有x，y，by，公共列表头名称gene_id；合并生成名为all_count的文件
all_count <- merge(merge(SRR3589960_count,SRR3589961_count,by="gene_id"),SRR3589962_count,by="gene_id")
#根据样本信息设置合并后矩阵列名依次为"gene_id","SRR60","SRR61","SRR62"
names(all_count)<-c("gene_id","SRR60","SRR61","SRR62")
head(all_count)
#gene_id SRR60 SRR61 SRR621
#__alignment_not_unique 1765232 1781423 11775682
#__ambiguous 32495 68314 749973
#__no_feature 430123 196950 1194454
#__not_aligned 46302661 31325937 385930145
#__too_low_aQual 4212734 3090224 34743276
#ENSG00000000003.15_7 0 0 0
all_count <- as.data.frame(all_count[-c(1:5),])查看all_count前5行（未合并时为count文件最后5行）发现并非所需的基因信息，可将其剔除；之后的分析主要依赖合并文件all_count进行。
2.3. Ensemble ID版本号的去除
即去除ENSG开头编号小数点后面的版本号，可调用stringr包的str_split函数进行数据分割。
##count文件名SRR3589960_count、列名gene_id
all_count$gene_id=unlist(str_split(all_count$gene_id,"[.]",simplify=T))[,1]2.4. Ensemble ID转换
借助getBM()函数进行ID转换之前，需通过数据库获取Ensemble_ID与Gene_ID之间的注释信息；选用数据库Ensembl中的人类基因注释信息hsapiens_gene_ensembl，由于注释信息繁多，这里只选择其中的ENS编号以及ENS编号对应的基因名称、染色体名称、起始与终止物理位置、转录本名称与长度。
library(biomaRt)
#listMarts()可列举其他可用数据库
human <- useMart('ensembl',dataset = "hsapiens_gene_ensembl")
##查看可注释信息，包括ensembl_gene id、基因名称、染色体名称、起始与终止物理位置、转录本名称与长度
listFilters(ensembl)
#提取ENS编号对应的注释信息（根据个需进行选择）
test <- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name",'start_position','end_position','ensembl_transcript_id','transcript_length'),mart = human)
head(test)
#对ensembl_gene_id升序排序
test <- test[order(test$ensembl_gene_id,decreasing = F),]
#去除重复的基因信息
test <- test[!duplicated(test$ensembl_gene_id),]
head(test)
#ensembl_gene_id hgnc_symbol chromosome_name start_position end_position ensembl_transcript_id transcript_length
#ENSG00000290166 19 1163451 1165333 ENST00000702095 489
#ENSG00000290165 X 119831616 119831759 ENST00000703415 144
#ENSG00000290164 X 119822906 119823043 ENST00000703414 138
#ENSG00000290163 X 119815164 119815307 ENST00000703413 144
#ENSG00000290162 X 119807433 119807581 ENST00000703412 149
#ENSG00000290149 7 37683843 37950412 ENST00000476620 555
此时的test注释文件仍有许多ENS编号没有对应的基因信息，可选择剔除hgnc_symbol列为空d的行
后续还需从这个总的注释文件中提取与我们count文件重叠的ENS编号及其对应的注释信息接下来将从以上过滤的注释信息中提取我们的count文件ENS编号(gene_id)所对应的注释信息，去除多余的无用信息。
思路：
1.提取all_count文件的gene_id(ENS编号),命名为id2.从注释文件test中提取id中的ENS编号及所在行的其他所有注释信息，命名为test_id3.对提取的注释信息进行过滤去重并排序4.添加列名5.合并注释信息与all_count矩阵6.导出含有基因名称等注释信息的count矩阵
代码如下
#取注释信息test与count矩阵的交集ensembl_id（共有的部分）
id=intersect(all_count$gene_id,test$ensembl_gene_id)
head(id)
#从注释文件test中提取共有id所在行的所有注释信息
test_id <- test[which(test$ensembl_gene_id%in%id),]
#ensembl_gene_id升序排序
test_id <- test_id[order(test_id$ensembl_gene_id,decreasing = F),]
##合并过滤后d的注释信息test_id与all_count，并注释列名
names(test_id)<-c("gene_id",'hgnc_symbol',"chromosome_name",'start_position','end_position','ensembl_transcript_id','transcript_length')
total_count <- merge(test_id,all_count,by='gene_id')
write.csv(total_count, "/Users/Zhanghp/Desktop/count_merge/mydata.csv")转换ID后的文件一览：
本部分完整代码：
##BiocManager::install("“biomaRt”")#biomart安装
##BiocManager::install("AnnotationDbi")
library(stringr)
library(biomaRt)
##载入count文件
SRR3589960_count <- read.table("/Users/Desktop/count_merge/SRR3589960.count",sep="\t", col.names = c("gene_id","control"))
SRR3589961_count <- read.table("/Users/Desktop/count_merge/SRR3589961.count",sep="\t", col.names = c("gene_id","control"))
SRR3589962_count <- read.table("/Users/Desktop/count_merge/SRR3589962.count",sep="\t", col.names = c("gene_id","control"))
##合并count文件
all_count <- merge(merge(SRR3589960_count,SRR3589961_count,by="gene_id"),SRR3589962_count,by="gene_id")
names(all_count)<-c("gene_id","SRR60","SRR61","SRR62")
all_count <- as.data.frame(all_count[-c(1:5),])
##去除gene_id中的版本号，即.之后的部分
all_count$gene_id=unlist(str_split(all_count$gene_id,"[.]",simplify=T))[,1]
listMarts()
human <- useMart('ensembl',dataset = "hsapiens_gene_ensembl")
listFilters(human)
test <- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name",'start_position','end_position','ensembl_transcript_id','transcript_length'),mart = human)
test <- test[order(test$ensembl_gene_id,test$transcript_length,decreasing = F),]
test <- test[!duplicated(test$ensembl_gene_id),]
head(test)
id=intersect(all_count$gene_id,test$ensembl_gene_id)
head(id)
test_id <- test[which(test$ensembl_gene_id%in%id),]
test_id <- test_id[order(test_id$ensembl_gene_id,decreasing = F),]
names(test_id)<-c("gene_id",'hgnc_symbol',"chromosome_name",'start_position','end_position','ensembl_transcript_id','transcript_length')
total_count <- merge(test_id,all_count,by='gene_id')
write.csv(total_count, "/Users/.../Desktop/mydata.csv")
支付宝扫一扫
微信扫一扫
.png)
.png)
