Извлечение количества прочтений, соответствующих каждому символу гена
Я провел количественную оценку экспрессии генов с помощью Salmon, которая дает мне транскрипты Ensembl, я преобразовал транскрипты Ensembl в символ гена, но для некоторых генов я использовал несколько транскриптов; Как я мог свернуть количество прочтений в гены, я попробовал tximport package, но обнаружил, что это слишком сложно, так как моя аннотация отличается.
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 в имя гена.
https://www.dropbox.com/s/m1iybfbu2i4bb39/Converting_transcript_id_to_gene_id.txt?dl=0
Возможная стратегия состоит в том, чтобы сгруппировать экспрессию генов по символам генов, а затем усреднить значения экспрессии по различным идентификаторам транскриптов Ensembl для одного и того же символа гена. Если бы вы предоставили репрезентативные данные выборки (ваши текущие данные выборки не содержат значений экспрессии из разных идентификаторов Ensembl ID с одинаковым символом гена), другим было бы намного проще помочь.
Пожалуйста, не используйте ссылки (кстати, ссылки не работают). Вместо этого создайте свои базовые данные и ожидаемый результат в своем запросе/вопросе. Достаточно пары строк. Чтобы люди, которые готовы помочь, получили точное представление о вашей проблеме.
Ответы 1
Вы можете использовать пакет dplyr.
Создайте тестовую таблицу:
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
Надеюсь это поможет.
Редактировать:
Как указано в комментариях ОП, ожидаемый результат поможет. Извините, может быть, я неправильно понял слово «крах» в этом контексте. Моя интерпретация заключается в суммировании прочтений на имя гена.
Редактировать2:
Как уже упоминалось в моем комментарии, постарайтесь запретить давать ссылки. Ссылки могут быть неработающими и т. д. Подробные инструкции о том, как написать хороший пост, смотрите: здесь.
Однако на основе ваших реальных данных сделайте следующее:
Загрузите данные:
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")
Просто объедините данные по этому идентификатору транскрипта:
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 читает больше всего.
Спасибо за то, что вы так полезны
Не могли бы вы добавить то, что вы ожидаете в качестве вывода, и какой код уже дал сбой?