feature R语言进行表达定量

最新推荐文章于 2024-04-17 21:20:29 发布

原创最新推荐文章于 2024-04-17 21:20:29 发布 · 496 阅读

8 ·

CC 4.0 BY-SA版权

文章标签：

#r语言 #数据库 #开发语言

文章介绍了如何使用R中的各种库（如Rsubread,limma,edgeR）处理bam和gtf/gff文件，进行转录本计数（featureCounts），并计算FPKM和TPM值，最后将统计结果和表达数据保存为表格。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

准备bam文件、gtf/gff文件

library(Rsubread)
library(limma)
library(edgeR)
library(tidyverse)
library(readr)
sample_name <- read_csv("sample.txt", col_names = FALSE) %>% pull(var = 1) %>% .[1:4] %>% str_sub(1,11)

for (i in 1:length(sample_name)) {
  
  sample = sample_name[i]
  bamFile = str_c("~/Program/tpm/",sample,".bam",sep = "")
  gtfFile = "~/Program/tpm/flax_trans.gtf"
  nthreads <- 16
  outFilePref <- "~/Program/tpm/"
  
  outStatsFilePath  <- str_c(outFilePref,sample, '.log',  sep = '')
  
  outCountsFilePath <- str_c(outFilePref,sample, '.count',  sep = '')
  
  fCountsList = featureCounts(bamFile,
                              annot.ext=gtfFile,
                              isGTFAnnotationFile=TRUE,
                              nthreads=nthreads,
                              isPairedEnd=F)
  dgeList = DGEList(counts=fCountsList$counts, genes=fCountsList$annotation)
  fpkm = rpkm(dgeList, dgeList$genes$Length)
  tpm = exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
  
  write.table(fCountsList$stat, outStatsFilePath, sep="\t", col.names=FALSE, row.names=FALSE, quote=FALSE)
  
  featureCounts = cbind(fCountsList$annotation[,1], fCountsList$counts, fpkm, tpm)
  colnames(featureCounts) = c('gene_id', 'counts', 'fpkm','tpm')
  write.table(featureCounts, outCountsFilePath, sep="\t", col.names=TRUE, row.names=FALSE, quote=FALSE)
  
}