如何绘制跨基因组坐标的 log2 倍数变化（使用 Deseq2 输出 csv）

Question

我有来自细菌基因组的 RNA-seq 数据（针对 2 种不同的处理重复 3 次），并使用 DeSeq2 计算了基因的 log2fc (padj < 0.05)。这会生成一个 csv 文件，其中包括（但不限于）基因名称和 log2fcexample of output.

更新：基因组已经发表并注释了，所以我有每个基因对应的基因组坐标。也许就像合并这些信息一样简单。但并不是所有的基因都是差异表达的，所以就比较复杂了...

但是，我想记录 2 RNA 变化（y 轴）与基因组坐标（x 轴）。但是我在互联网上搜索没有成功。有谁知道一个相对简单的方法来做到这一点？我很高兴使用 R/python...我已经包含了一篇论文中的示例以说明我所追求的...example of what I'm after

也许这就是为什么没有人谈论它的原因。但是在我附上的图片中，他们没有讨论他们是如何绘制的。

提前致谢！！

Answer 1

让我们用一个例子来获得 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))

我们需要获取一个注释，可以是床文件，你使用来自 rtracklayer 的相同功能import()：

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>

您需要将 deseq2 结果的行名与列匹配，在本例中为 gene_id，如果您的 table 是符号，请将其与适当的匹配：

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

这里有不止一条染色体，所以我们只做 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"))

如何绘制跨基因组坐标的 log2 倍数变化（使用 Deseq2 输出 csv）

How can I plot log2 fold-change across genome coordinates (using Deseq2 output csv)

python

r

rna-seq