Saya memiliki data RNA-seq (3 ulangan untuk 2 perlakuan berbeda) dari genom bakteri dan telah menggunakan DeSeq2 untuk menghitung log2fc untuk gen (padj <0,05). Ini menghasilkan file csv yang menyertakan (tetapi tidak terbatas pada) nama gen dan log2fccontoh keluaran.

Pembaruan: Genom diterbitkan dan diberi anotasi, jadi saya memiliki koordinat genom yang sesuai untuk setiap gen. Mungkin sesederhana menggabungkan informasi ini. Tetapi tidak semua gen diekspresikan secara berbeda, sehingga menjadi lebih rumit...

Namun, saya ingin mencatat perubahan RNA log2 (sumbu y) vs koordinat genom (sumbu x). Tapi saya menjelajahi internet tanpa hasil. Adakah yang tahu cara yang relatif sederhana untuk melakukan ini? Saya senang menggunakan R/python... Saya telah menyertakan contoh dari makalah untuk apa yang saya cari... contoh apa yang saya cari

Mungkin ini sangat sederhana sehingga tidak ada yang membicarakannya. tetapi dalam gambar yang saya lampirkan, mereka tidak membahas bagaimana mereka merencanakannya.

Terima kasih sebelumnya!!

0
neon8731 5 Juli 2020, 11:25

1 menjawab

Jawaban Terbaik

Mari kita gunakan contoh untuk mendapatkan hasil DESeq2:

library(DESeq2)
library(zebrafishRNASeq)
data(zfGenes)

head(zfGenes)
                   Ctl1 Ctl3 Ctl5 Trt9 Trt11 Trt13
ENSDARG00000000001  304  129  339  102    16   617
ENSDARG00000000002  605  637  406   82   230  1245
ENSDARG00000000018  391  235  217  554   451   565
ENSDARG00000000019 2979 4729 7002 7309  9395  3349

dds = DESeqDataSetFromMatrix(zfGenes[rowMeans(zfGenes)>30,],
data.frame(grp=gsub("[0-9]*","",colnames(zfGenes))),~grp)
dds = DESeq(dds)

res = data.frame(results(dds))

Kita perlu mendapatkan anotasi, bisa berupa file tempat tidur, Anda menggunakan fungsi yang sama import() dari rtracklayer:

library(rtracklayer)
gtf = import("ftp://ftp.ensembl.org/pub/release-100/gtf/danio_rerio/Danio_rerio.GRCz11.100.gtf.gz")

genes = gtf[gtf$type=="gene"]

head(genes)
GRanges object with 6 ranges and 20 metadata columns:
      seqnames      ranges strand |   source     type     score     phase
         <Rle>   <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
  [1]        4 17308-18211      - |  ensembl     gene      <NA>      <NA>
  [2]        4 31259-45642      + |  ensembl     gene      <NA>      <NA>
  [3]        4 38344-45396      + |  ensembl     gene      <NA>      <NA>
  [4]        4 48519-53370      - |  ensembl     gene      <NA>      <NA>
  [5]        4 57034-64703      - |  ensembl     gene      <NA>      <NA>
  [6]        4 69619-72011      + |  ensembl     gene      <NA>      <NA>

Anda harus mencocokkan nama baris hasil deseq2 dengan kolom, dalam hal ini gene_id, jika tabel Anda dalam simbol, cocokkan dengan yang sesuai:

annotation = as.data.frame(genes)[match(rownames(res),genes$gene_id),]
head(as.data.frame(genes)[match(rownames(res),genes$gene_id),])
      seqnames    start      end width strand         source type score phase
12565        9 34112067 34121839  9773      - ensembl_havana gene    NA    NA
12564        9 34089156 34113209 24054      + ensembl_havana gene    NA    NA
484          4 15081385 15103696 22312      - ensembl_havana gene    NA    NA
480          4 15011341 15059876 48536      + ensembl_havana gene    NA    NA
21778       12 33484458 33537126 52669      + ensembl_havana gene    NA    NA

Ada lebih dari satu kromosom di sini jadi mari kita lakukan kromosom 25:

res = cbind(res,annotation[,1:5])
res = res[which(res$seqnames == 25),]

plot(res$start,res$log2FoldChange,xaxt="n",xlab="Pos on chr25",ylab="Log2FC")
axis(side=1,5e6*1:8,paste(5 * 1:8,"Mb"))

enter image description here

0
StupidWolf 5 Juli 2020, 14:32