把基因画在染色体上

2020-06-11 本文已影响0人果蝇饲养员的生信笔记

之前有一次，我导说“能不能把这些基因按顺序画在果蝇的染色体上”，因此，我寻找了一些实现这个需求的方法。经过一番奇奇怪怪的操作之后，最终用ggplot2画了。。。其它方法简单地记录在这里（反正也是学了的，万一以后有用呢）

生物实验室的日常-表情包

一、画基因和染色体的R包

RIdeogram：画染色体很好看
gggenes：画基因的结构
ggbio：可以把基因画在染色体上
Gviz：把基因组坐标和其它数据用相同的坐标放在一起展示

用处不大，随便看看

二、安装R包

options("repos" = c(CRAN="https://mirror.lzu.edu.cn/CRAN/")) 
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor")
install.packages("RIdeogram")
install.packages("gggenes")
BiocManager::install("ggbio")
BiocManager::install("Gviz")
library(RIdeogram)
library(RColorBrewer)
library(ggbio)
library(gggenes)
library(Gviz)
library(ggplot2)

三、RIdeogram

（1）每条染色体的长度

需要至少有Chr/ Start/ End三列，如果还有着丝粒位点可以加上CE_start/ CE_end两列。可以从UCSC查到所需物种的每条染色体长度，这里的染色体编号要和后面的一致。

> DmChromosome <- read.table("DmChromosome.txt", header=TRUE, sep="\t", row.names=1, stringsAsFactors=F) 
> DmChromosome 
  Chr Start      End
1  2L     0 23513712
2  2R     0 25286936
3  3L     0 28110227
4  3R     0 32079331
5   4     0  1348131
6   X     0 23542271
7   Y     0  3667352

（2）基因密度信息

可以通过gffex函数从gff文件中提取热图信息（基因密度）。多次尝试选择一个合适的window大小，这里的feature也可以是CDS、exon、ncRNA等。

> Dmdensity <- GFFex(input="Drosophila_melanogaster.BDGP6.28.99.chr.gff3.gz",
                   karyotype="DmChromosome.txt",
                   feature="gene",window=200000)
> head(Dmdensity)
  Chr   Start     End Value
1  2L       1  200000    25
2  2L  200001  400000    29
3  2L  400001  600000    31
4  2L  600001  800000    14
5  2L  800001 1000000    23

（3）marker型绘图

有marker/line/polygon/heatmap四种绘图方法。marker型需要Type/ Shape/ Chr/ Start/ End/ color。它仅对type一栏绘制图例，marker型只适合展示少量基因。

> head(gene_list) 
  Type Shape Chr    Start      End  color   name
1 mRNA   box  2L 18683720 18688800 FFFFB3 gene-1
2 mRNA   box  2L 18684984 18688800 FFFFB3 gene-1
3 mRNA   box  3L  7123213  7125533 FFFFB3 gene-3
4 mRNA   box  2L  3462218  3466112 FFFFB3 gene-4
5 mRNA   box  2L  3462218  3466112 FFFFB3 gene-4
6 mRNA   box   X  5873840  5876853 FB8072 gene-2
> ideogram(karyotype=DmChromosome, #染色体信息
+          overlaid=Dmdensity, #基因密度
+          label=gene_list, #目标基因标签
+          label_type="marker",width=135)
> convertSVG("chromosome.svg",device="png") #非常耗时

确实是好看呢

（4）line/polygon/heatmap型绘图

需要制作一个类似于基因密度信息的表，包含Chr/ Start/ End/ Value_1/ Color_1/ Value_2/ Color_2。

> head(gene_list_v)
  Chr   Start     End      Value  Color
1  2L       1  200000  1966.7804 fc8d62
2  2L  200001  400000  4279.0219 fc8d62
3  2L  400001  600000  3022.2958 fc8d62
4  2L  600001  800000   182.4522 fc8d62
5  2L  800001 1000000 14683.3203 fc8d62
6  2L 1000001 1200000  1023.5985 fc8d62
> ideogram(karyotype=DmChromosome, #染色体信息
+          overlaid=Dmdensity, #基因密度
+          label=gene_list_v, #目标基因标签
+          label_type="line",width=135)
> convertSVG("chromosome.svg",device="png") #可以手动修改文件名

line

（5）heatmap型绘图

需要制作一个类似于基因密度信息的表，包含Chr/ Start/ End/ Value。

ideogram(karyotype=DmChromosome,
         overlaid=Dmdensity,
         label=gene_list_v,
         label_type="heatmap",width=135,
         colorset1 = c("#f7f7f7", "#2c7fb8"),
         colorset2 = c("#f7f7f7", "#e34a33"))
convertSVG("chromosome.svg",device="png")

heatmap

（6）数据集

如果是人类的染色体，可以通过直接获得核型信息（不需要这里的前两步）。

data(human_karyotype,package="RIdeogram")
data(gene_density,package="RIdeogram") 
head(human_karyotype) #24行 Chr/ Start/ End/ CE_start/ CE_end
head(gene_density) #3102行 Chr/ Start/ End/ Value

四、gggenes

（1）示例数据

> head(example_genes); dim(example_genes)
  molecule gene  start    end  strand direction
1  Genome5 genA 405113 407035 forward         1
2  Genome5 genB 407035 407916 forward         1
3  Genome5 genC 407927 408394 forward         1
4  Genome5 genD 408387 408737 reverse        -1
5  Genome5 genE 408751 409830 forward         1
6  Genome5 genF 409836 410315 forward         1
[1] 72  6

（2）绘图

> p1 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule,fill=gene)) +
+   geom_gene_arrow() +
+   facet_wrap(~ molecule,scales="free",ncol=1) +
+   scale_fill_brewer(palette="Set3") +
+   theme_genes()
> p1

（3）对齐某一基因

> dummies <- make_alignment_dummies(example_genes,
+                                   aes(xmin=start,xmax=end,y=molecule,id=gene),
+                                   on="genE")
> p1 + geom_blank(data=dummies)

（4）写上基因名称

> p2 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule,fill=gene)) +
+   geom_gene_arrow(arrowhead_height=unit(3,"mm"),arrowhead_width=unit(1,"mm")) +
+   facet_wrap(~ molecule,scales="free",ncol=1) +
+   scale_fill_brewer(palette="Set3") +
+   theme_genes()
> p2
> p2 + geom_gene_label(aes(label=gene),align="left")

（5）区分正负链

> p3 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule,
+                                fill=gene,forward=direction)) +
+   geom_gene_arrow() +
+   facet_wrap(~ molecule,scales="free",ncol=1) +
+   scale_fill_brewer(palette="Set3") +
+   theme_genes()
> p3

（6）基因亚结构

> head(example_subgenes);dim(example_subgenes)
  molecule gene  start    end  strand subgene   from     to
1  Genome5 genA 405113 407035 forward  genA-1 405774 406538
2  Genome5 genB 407035 407916 forward  genB-1 407458 407897
3  Genome5 genC 407927 408394 forward  genC-1 407942 408158
4  Genome5 genC 407927 408394 forward  genC-2 408186 408209
5  Genome5 genC 407927 408394 forward  genC-3 408233 408257
6  Genome5 genF 409836 410315 forward  genF-1 409938 410016
[1] 143   8
> p4 <- ggplot(example_genes,aes(xmin=start,xmax=end,y=molecule)) +
+   facet_wrap(~ molecule, scales="free",ncol = 1) +
+   geom_gene_arrow(fill="white") +
+   geom_subgene_arrow(data=example_subgenes,
+                      aes(xmin=start,xmax=end,y=molecule,fill=gene,
+                          xsubmin=from,xsubmax=to),color="black",alpha=.7) +
+   theme_genes()
> p4

五、ggbio

只用到了一个功能，把一些基因（比如上调或下调的基因）画在染色体上。

（1）构建GRanges对象

> mygr #GRanges object
GRanges object with 127 ranges and 4 metadata columns:
        seqnames          ranges strand | gene_name        logFC      pvalue regulation
           <Rle>       <IRanges>  <Rle> |  <factor>    <numeric>   <numeric>   <factor>
    [1]     chrX 3720007-3788818      + |       Tlk -6.982227392 0.394290308       down
    [2]     chrX 3858940-3862697      - |      roX1  11.48670351  0.00280582         up
    [3]     chrX 3932156-3936469      - |    CG2938 -4.634468751 0.373950298       down

（2）绘图

> autoplot(mygr,layout="karyogram",aes(color=regulation,fill=regulation))
> mygr$levels <- as.numeric(factor(mytx$regulation))
> autoplot(mygr,layout="karyogram",aes(color=regulation,fill=regulation,
+                                      ymin=(levels-1)*10/2,ymax=levels*10/2))

六、Gviz

Gviz能展示的内容非常多。特点是把基因组坐标和其它数据（如logFC）用相同的坐标放在一起展示。每一部分内容作为一个track，可以有GenomeAxisTrack、IdeogramTrack、DataTrack、AnnotationTrack、GeneRegionTrack、 SequenceTrack、 AlignmentsTrack等tracks。
（感觉过于复杂，我又不太能用到，所以只是粗略地看了看。官网的文档很详细）

（1）构造GRanges对象

> mygrsub[1:3] 
GRanges object with 3 ranges and 6 metadata columns:
      seqnames          ranges strand |   symbol          CF          CM        logFC
         <Rle>       <IRanges>  <Rle> | <factor>   <numeric>   <numeric>    <numeric>
  [1]    chr2L 6911090-6914138      - |     Wee1 11.99432931 1.329478214 -3.173420446
  [2]    chr2L 6911343-6914138      - |     Wee1 4.959662983  68.0307017  3.777871974
  [3]    chr2L 6914301-6920210      - |      x16 3.455682714 2.584117733  -0.41929896
           pvalue       group
        <numeric> <character>
  [1] 0.066684944        Wee1
  [2] 0.000146948        Wee1
  [3] 0.644895204         x16
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

（2）绘图

> g <- GenomeAxisTrack()
> a <- AnnotationTrack(mygrsub,name="gene ranges",feature=updown)
> d <- DataTrack(mygrsub, data="logFC", baseline=0,
+                type="h", name="log2 fold change", strand="+") #h,histogram lines
> displayPars(g) <- list(col="darkgray",fontcolor="lightblue4")
> displayPars(a) <- list(background.panel="#FFFEDB",background.title="brown",col=NULL)
> displayPars(d) <- list(background.panel="azure",background.title="dodgerblue4",col=NULL)
> plotTracks(list(g,d,a),groupAnnotation="group",
+            down="#40B3D9", up="#FE8C8C")

最上面的数轴代表GenomeAxisTrack，log2FC就放在了DataTrack，基因区间放在了AnnotationTrack。其中DataTrack就是想要展示的数据了。还可以有很多种绘图方法，像dot and lines plot、 boxplot、 heatmap等。

文档很详细

把基因画在染色体上

一、画基因和染色体的R包

二、安装R包

三、RIdeogram

（1）每条染色体的长度

（2）基因密度信息

（3）marker型绘图

（4）line/polygon/heatmap型绘图

（5）heatmap型绘图

（6）数据集

四、gggenes

（1）示例数据

（2）绘图

（3）对齐某一基因

（4）写上基因名称

（5）区分正负链

（6）基因亚结构

五、ggbio

（1）构建GRanges对象

（2）绘图

六、Gviz

（1）构造GRanges对象

（2）绘图

猜你喜欢

热点阅读