提取与每个基因符号匹配的读数

Question

我已经通过 Salmon 量化了基因表达，它为我提供了 Ensembl 转录本，我将 Ensembl 转录本转换为基因符号，但对于某些基因，我使用了多个转录本；我如何将读取计数折叠到基因中，我尝试了 tximport 包，但我发现这太难了，因为我的注释不同。

Name                NumReads
ENST00000355520.5   407.186
ENST00000566753.1   268.879
ENST00000481617.2   242.25
ENST00000538183.2   226.576

ensembltranscript_id gene_name
ENST00000482226.2   FCGR2C
ENST00000508651.1   FCGR2C
ENST00000571914.1   TSPAN10
ENST00000571707.1   TSPAN10
ENST00000534817.1   OVCH2
ENST00000445557.1   OR52E1
ENST00000575319.1   CYP2D7
ENST00000576465.1   CYP2D7

EDITED

这是 Salmon 读取计数的输出

https://www.dropbox.com/s/7bkril0v6sw7v9z/Salmon_output.txt?dl=0

这是我将 Salmon 输出中的转录本 ID 转换为基因名称的时候

https://www.dropbox.com/s/m1iybfbu2i4bb39/Converting_transcript_id_to_gene_id.txt?dl=0

Answer 1

您可以使用包 dplyr。

创建测试table:

names = c("ensembltranscript_id",  "gene_name", "NumReads")
transcripts = c("ENST00000482226.2", "ENST00000508651.1", "ENST00000571914.1", "ENST00000571707.1", "ENST00000534817.1")
gene_names = c("FCGR2C", "FCGR2C", "TSPAN10", "TSPAN10", "OVCH2")
reads = c(205.56, 456.21, 123.3, 52.6, 268.45)

data = data.frame(transcripts, gene_names, reads)
names(data) = names

进行计算：

result = data %>% 
  group_by(gene_name) %>% 
  summarise(sum(NumReads)) %>%
  mutate_if(is.numeric, format, 2)

打印结果：

# A tibble: 3 x 2
  gene_name `sum(NumReads)`
  <fct>     <chr>          
1 FCGR2C    661.77         
2 OVCH2     268.45         
3 TSPAN10   175.90

希望这对您有所帮助。

编辑：

正如 OP 的评论中所述，预期的输出会有所帮助。抱歉，也许我在这方面误解了 'collapse'。我的解释是将每个基因名称的读数相加。

编辑2：

正如我评论中提到的，尽量避免提供链接。链接可能会损坏等。有关如何编写好的 post 的完整说明，请参阅：here。

但是，根据您的真实数据执行以下操作：

加载数据：

salmon_reads = read.table(file = "/path/to/Salmon_output.txt", header = T, sep = "\t")
genes = read.table(file = "/path/to/Converting_transcript_id_to_gene_id.txt", header = T, sep = "\t")

只需将数据合并到 transcript-id:

merged_data = merge(x = salmon_reads, y = genes, by.x = colnames(salmon_reads)[1], by.y = colnames(genes)[1], all = T)

进行减少读数的计算和排序：

result = merged_data %>% 
  group_by(external_gene_name) %>% 
  summarise(sum(NumReads)) %>%
  mutate_if(is.numeric, format, 2)

result$`sum(NumReads)` = as.numeric(result$`sum(NumReads)`)
result = result[order(result$`sum(NumReads)`, decreasing = T),]

你没有提到如何处理 NA。在这种情况下，汇总所有 NA 基因名称的读数。这就是为什么 NA 的阅读量最多的原因。

提取与每个基因符号匹配的读数

Extracting read counts matched to each gene symbol

loops

r

rna-seq