Цель
Я пытаюсь получить распределение размеров экзонов и интронов у трехиглой колюшки ( Gasterosteus aculeatus ) .
Данные загружены
Я скачал некоторые данные из Ensembl. Точнее я туда зашел , в "атрибутах" выбрал "структура" и выбрал
- Ensembl Gene ID
- Exon Chr Start (bp)
- Exon Chr End (bp)
В этих данных есть две проблемы:
Первый пункт я предположил из-за реверсии, поэтому я просто поменял местами начальную и конечную позиции, когда такая инверсия была. Второй пункт более проблематичен, поскольку я не знаю, что такое совпадение может означать на самом деле.
Не могли бы вы помочь найти распределение размеров экзонов и интронов у трехиглой колюшки?
Что я сделал с приведенными выше данными
Подготовка
# read data
d = read.table(file.choose(), header=TRUE)
names(d)= c("ID", "start", "end")
# Reverse exons when needed
toReverse = which(d$start > d$end)
s = d$start[toReverse]
d$start[toReverse] = d$end[toReverse]
d$end[toReverse] = s
Размеры экзонов -> Выглядит хорошо
# Exon sizes
ExonSizes = d$end - d$start
hist(ExonSizes) # Cool, looks good!
Размеры интронов -> Выглядит плохо
# Intron sizes
# This is a little slow. One might have better to switch to C++.
# The idea is to look at the distance between the end of an exon and the start of the next exon within each gene.
d$ID = paste(d$ID)
IntronSizes = c()
uniquegenes = unique(d$ID)
nbgenes = length(uniquegenes)
i=0
exon = 0
for (gene in uniquegenes)
{
i=i+1
cat(paste0(i," / ",nbgenes,"\n"))
while (TRUE)
{
exon=exon+1
if (d$ID[exon] == gene) {
IntronSizes = append(IntronSizes, d[exon+1,]$start - d[exon,]$end)
} else
{
break
}
}
}
hist(IntronSizes) # There are negative values. I don't know how to interpret them.
hist(log(abs(IntronSizes))) # That looks good but I doubt it makes much sense!
В командной строке:
wget ftp://ftp.ensembl.org/pub/release-84/gtf/gasterosteus_aculeatus/Gasterosteus_aculeatus.BROADS1.84.gtf.gz
gunzip Gasterosteus_aculeatus.BROADS1.84.gtf.gz
Затем в Р:
library(GenomicFeatures)
txdb = makeTxDbFromGFF("Gasterosteus_aculeatus.BROADS1.84.gtf")
hist(log10(width(unique(exons(txdb))))) # exons
hist(log10(width(unique(unlist(intronsByTranscript(txdb)))))) # introns
Обратите внимание, что некоторые из аннотированных экзонов (и интронов между ними) вряд ли будут правильными. Например, есть экзоны и интроны длиной 1 основание. Я не думаю, что кто-то на самом деле верит, что это правильно, но в остальном распределение, вероятно, приблизительно правильное.
Редактировать : для чего это стоит, в файле GTF нет экзонов отрицательного размера. Я предполагаю, что перекрывающиеся экзоны происходят из разных транскриптов (или генов на противоположных цепях).
Edit2 : Если вы хотите получить свои интроны из «модели союзного гена», используйте что-то вроде этого, reduce(exonsBy(txdb, by='gene'))
а затем lapply
gaps()
к этому. Результаты будут правильными, и процесс, вероятно, займет меньше времени, чем то, что вы пробовали.
GenomicFeatures
это пакет биокондуктора и что в MAC OSX его wget
необходимо заменить наcurl
фон Мизес