Я ищу способ извлечь attributes
Note
и ID
информацию из файла GFF. Я загрузил файл GFF, используя
library(ape)
GFF<-ape::read.gff(file = "gene_models.gff")
. Я хочу сделать что-то подобное, чтобы сопоставить набор данных (созданный из DESeq2
) со столбцом атрибутов Note и ID:
Data <- GFF$attributes(ID)[match(DESeq2$seqid, GFF$attributes(Note))]
DESeq2
было выполнено:
salmon quant -i salmon_index -l A -1 input_R1_paired.fastq -2 input_R2_paired.fastq -o transcripts
library(GenomicFeatures)
library(tximport)
library(DESeq2)
library(AnnotationDbi)
dir <- "/path"
list.files(dir)
samples <- read.csv2("path",sep =";", header = TRUE)
head(samples)
my.files <- list.files(list.dirs(path = ".", full.names = TRUE, recursive = FALSE), pattern = "quant.sf", full.names = TRUE)
my.files
names(my.files) <- sapply(strsplit(my.files, split="\\/"), function(x)x[2])
head(my.files, 6)
all(file.exists(my.files))
TxDb <- makeTxDbFromGFF(file = "gene_models.gff")
k <- keys(TxDb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(TxDb, k, "GENEID", "TXNAME")
head(tx2gene)
txi <- tximport(files = my.files, type = "salmon", tx2gene = tx2gene)
head(txi$counts)
rownames(samples)<-samples$Run
colnames(txi$counts)
rownames(samples) <- colnames(txi$counts)
samples1 <- data.frame(Group= factor(rep(c("High", "Low"), each = 3)))
rownames(samples1) <- colnames(txi$counts)
ddsTxi <- DESeqDataSetFromTximport(txi, colData = samples1, design = ~ Group)
ddsTxi
ddsTxi <- DESeq(ddsTxi)
res <- results(ddsTxi)
head(results(ddsTxi, tidy=TRUE))
summary(res)
res <- res[order(res$padj),]
head(res)
write.csv(as.data.frame(res), file="res.csv")
И выходной файл выглядит так:
seqid baseMean log2FoldChange lfcSE stat pvalue padj
<character> <numeric> <numeric> <numeric> <numeric> <numeric <numeric>
1 G_06g024400.1 1266.29960801482 -3.27190177919055 0.493708714717401 -6.62719065241412 3.42135630989487e-11 6.20155044731545e-07
И GFF
выглядит так:
seqid source type start end score strand phase attributes
1 G1.3ch00 maker gene 21882 22007 NA - <NA> ID=G_06g024400.1;Name=G_06g024400.1;Note=Similar to 06g024400.1.1: LOW QUALITY:protein;ref_id=06g024400.1.1;
Спасибо за заранее!