背景
RNA-seq read count转换成TPM。
应用场景
只能拿到RNA-seq数据的read count，想转换成TPM。
注：TPM的计算方式最好是直接从原始的FASTQ文件进行获取，例如Samlon, Sailfish 或 kallisto, 或者先用FASTQ比对得到BAM，然后用RESM估计TPM.
count、fpkm和tpm对比介绍，参考:
https://bioinfoer.com/forum.php?mod=viewthread&tid=551&extra=page%3D1
下载TCGA RNA-seq的read count数据
如果你自己的RNA-seq数据已经保存成 easy_input.txt的格式，就跳过这步，直接进入“输入数据”
此处选择旧版本tcga下载得到的4个样本的conut数据。
新版本count和tpm下载方式参考：
https://bioinfoer.com/forum.php?mod=viewthread&tid=548&extra=page%3D1
输入数据
第一列是基因ID，第一行是sample ID。
一个单元格内是一个基因在一个sample中的read count。
#expMatrix <- read.table("easy_input.txt",
# row.names = 1, header = TRUE, sep="\t",as.is = T)
str(expMatrix)
#查看前三个基因的read count
expMatrix[1:3,]
计算基因长度
这一步有两种思路:
- 方法1：简单粗暴的计算基因在染色体的起始和结束之差
- （推荐）方法2：比较麻烦的计算每个基因的最长转录本(外显子之和)
方法1
如果只是想粗略了解一下表达情况，可以简单把基因在染色体上的起始位置和结束位置之差用作标准化的长度。
#source("https://bioconductor.org/biocLite.R")
#biocLite("biomaRt")
library(biomaRt)
#查看基因组参数
mart = useMart('ensembl')
listDatasets(mart)
#你需要哪个基因组，就复制它在dataset列里的词，放在下面这行的`dataset = `参数里
#此处以人类为例，植物参考注一
bmart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "hsapiens_gene_ensembl",
host = "www.ensembl.org")
# 从输入数据里提取基因名
feature_ids <- rownames(expMatrix)
attributes = c(
"ensembl_gene_id",
#"hgnc_symbol",
"chromosome_name",
"start_position",
"end_position"
)
filters = "ensembl_gene_id"
feature_info <- biomaRt::getBM(attributes = attributes,
filters = filters,
values = feature_ids, mart = bmart)
mm <- match(feature_ids, feature_info[[filters]])
feature_info_full <- feature_info[mm, ]
rownames(feature_info_full) <- feature_ids
# 计算基因的有效长度
eff_length <- abs(feature_info_full$end_position - feature_info_full$start_position)
names(eff_length) <- feature_info_full$ensembl_gene_id
这里的 eff_length就是每个基因的有效长度
方法2
这一步原本打算用Bioconductor的 ord.Hs.eg.db,但是实际发现里面的ENSEMBL ID明显比输入的TCGA数据集里的ID少。所以决定从GTF文件提取：
先写一个函数，实现从下载GTF文件到计算基因长度
get_eff_len <- function(file=NULL,
feature = "exon",
gene_id = 'gene_id',
transcript_id = "transcript_id",...){
# check if file parameter is provided
if ( is.null(file)){
stop(" no file path or url provided")
}
file_split <- unlist(strsplit(file,'/'))
file_name <- file.path(Sys.getenv('R_USER'),
file_split[length(file_split)])
# check file is url or local path
if (grepl(pattern = '(ftp|http|https)', file)){
file_name <- file_name
if ( ! file.exists(file_name)){
print(paste0("file will save in ", file_name))
download.file(url = file,
destfile = ,file_name,
method = "internal")
}
} else{
if ( ! file.exists(file)){
stop("file is not exists")
}
file_name <- file
}
# get the file type from the file name
# and build connection with different read function
file_type <- sub('.*?\\.(.*?)$','\\1',file_name)
if (file_type == 'gtf'){
gtf_in <- file(file_name, open='rt')
} else if(file_type == 'gz'){
file_name_tmp <- gsub(sprintf("[.]%s$",'gz'),"",file_name)
if ( ! file.exists(file_name_tmp))
R.utils::gunzip(file_name, remove = FALSE)
gtf_in <- file(file_name_tmp, open='rt')
} else{
err_message <- paste0("unsupport file type: ", file_name)
stop(err_message)
}
gtf_lines <- readLines(gtf_in, warn=FALSE)
close(gtf_in)
# save the effitve length of gene
feature_pattern <- paste0("\\t",feature,"\\t")
mt_rows <- sum(grepl(feature_pattern, gtf_lines))
# check the feature
if ( mt_rows <= 1){
features <- do.call(rbind, strsplit(gtf_in, '\\t'))[,3]
print("feature is unfoundable, the top 5 feature is:")
print(as.data.frame(table(feature))[1:5,])
stop(paste0("set your feature and run again"))
}
mt <- matrix(nrow = mt_rows, ncol= 1)
factor_names <- vector(mode = "character",length = mt_rows)
# save the gene id of the previous line
# iteration
line_num <- 1
for (line in gtf_lines){
if ( grepl(pattern = "^#", x = line)) next
records <- unlist(strsplit(x = line, '\\t'))
if ( grepl(feature, records[3]) ){
gene_pattern <- paste0(".*?", gene_id, " \"(.*?)\";.*")
gene_id_ <- sub(gene_pattern,"\\1",records[9])
tx_pattern <- paste0(".*?", transcript_id, " \"(.*?)\";.*")
tx_id_ <- sub(tx_pattern, "\\1",records[9])
gene_start <- records[4]
gene_end <- records[5]
# add 1 for gtf coordinate is 1-based
gene_len <- abs(as.numeric(gene_end) - as.numeric(gene_start) + 1)
mt[line_num,] <- gene_len
factor_names[line_num] <- paste0(gene_id_,"_",tx_id_)
line_num <- line_num + 1
}
}
by_transcript <- as.factor(factor_names)
tx_length <- sapply(split(mt[,1], by_transcript), sum, na.rm=TRUE)
by_gene <- as.factor(do.call(rbind,
strsplit(names(tx_length), '_'))[,1])
gene_length <- sapply(split(tx_length, by_gene), max)
return(gene_length)
}
目前只测试了拟南芥和人类。其他物种报错或者没有返回预期结果，请将物种的URL和你遇到的问题，请回复留言。
注：函数不支持需要登陆才能下载的GTF的URL
下面调用函数下载GTF文件，提取存放有效转录本长度的向量
把你想要的gtf文件的地址写在 url <- 的后面
案例一：
此处以TCGA的read count作为输入，用跟TCGA一致的注释文件提取外显子长度。
到这里查询TCGA用的是哪个版本的注释文件：https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files
旧版本有注释用的是v22：GDC.h38 GENCODE v22 GTF (used in RNA-Seq alignment and by HTSeq)，但是新版本已经没有，目前是star软件比对的。这里先用v22演示。刚好count文件也是旧版本的tcga下载的。
去刚才的页面下载v22版本gtf文件，解压缩到当前文件夹。
#url <- "ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_22/gencode.v22.annotation.gtf.gz"
file<- "gencode.v22.annotation.gtf"
eff_length <- get_eff_len(file = file, #或者换成file = url，自动下载解压缩
feature = "exon", #注意，用CDS时不会统计非编码RNA
gene_id = 'gene_id',
transcript_id = "transcript_id")
eff_length1 <- data.frame(eff_length)
write.csv(eff_length1, "eff_length.csv", row.names = TRUE)
eff_length1前几行如下：
gtf解压后文件较大，后续资料不提供，但是提供了转换后的“eff_length.csv”，不想转换的，可以跳过直接使用。
案例二：（未测试）拟南芥的GTF文件(推荐Araport)，统计每个基因最长的转录本的CDS序列长度
file <- "Araport11_GFF3_genes_transposons.201606.gtf"
eff_length <- get_eff_len(file = file,
feature = "CDS", #只看蛋白编码基因
gene_id = 'gene_id',
transcript_id = "transcript_id")
发表于 昨天 15:18
|
查看: 6|
回复: 2
|Archiver|手机版|小黑屋|bioinfoer
( 萌ICP备20244422号 )
