对测序得到的reads进行计数，即基因表达的定量过程。根据reads和基因位置的overlap，以此来判断reads到底属于哪一个基因，同时对该reads总数进行计数，生成counts矩阵。
今日内容
-count对reads进行计数
语言完成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 input
-
featuresfilename Path to the file containing the features
-
-
optional 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 be
-
sorted either by position or by read name, and the
-
sorting order must be specified. Ignored for single-
-
end data.
-
输入文件的排序方式，默认按read名排序
-
--max-reads-in-buffer MAX_BUFFER_SIZE
-
When <alignment_file> is paired end sorted by
-
position, allow only so many reads to stay in memory
-
until the mates are found (raising this number will
-
use more memory). Has no effect for single end or
-
paired 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 strand
-
interpretation
-
-a MINAQUAL, --minaqual MINAQUAL
-
skip all reads with alignment quality lower than the
-
given minimum value (default 10）
-
剔除mapping quality值低于阈值的read
-
-t FEATURETYPE, --type FEATURETYPE
-
feature type (3rd column in GFF file) to be used, all
-
features of other type are ignored (default, suitable
-
for Ensembl GTF files: exon)
前一期利用SAMtools工具对比对产生的sam文件进行格式转换并排序，得到了按染色体位置进行排序的bam文件${id}.。
因此输入文件格式-f bam，文件名${id}.sor ;排序方式-r pos ;文献中提到reads计数时剔除了mapping quality < 30的低质量reads，即-a 30。
输入bam文件路径：./aligned/${id}.
参考基因组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}. \
-
./genome_gtf/gencode. >./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包可通过'("package name")'命令进行下载安装。但使用以上命令安装biomaRt及其依赖包AnnotationDbi时将报错，无法获取安装包信息，可使用BioMannager进行安装，亲测有效
-
('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, = by, all = FALSE, = all, = all,
-
sort = TRUE, suffixes = c(".x",".y"), = 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 <- (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_id 3.对提取的注释信息进行过滤去重并排序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')
-
(total_count, "/Users/Zhanghp/Desktop/count_merge/")
转换ID后的文件一览：
本部分完整代码：
-
##BiocManager::install("“biomaRt”")#biomart安装
-
##BiocManager::install("AnnotationDbi")
-
library(stringr)
-
library(biomaRt)
-
##载入count文件
-
SRR3589960_count <- ("/Users/Desktop/count_merge/",sep="\t", = c("gene_id","control"))
-
SRR3589961_count <- ("/Users/Desktop/count_merge/",sep="\t", = c("gene_id","control"))
-
SRR3589962_count <- ("/Users/Desktop/count_merge/",sep="\t", = 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 <- (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')
-
(total_count, "/Users/.../Desktop/")