把基因画在染色体上

之前有一次，我导说“能不能把这些基因按顺序画在果蝇的染色体上”，因此，我寻找了一些实现这个需求的方法。经过一番奇奇怪怪的操作之后，最终用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等。


文档很详细
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 195,585评论 5赞 462
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 82,283评论 2赞 373
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 142,760评论 0赞 324
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,461评论 1赞 266
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,280评论 4赞 357
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,268评论 1赞 273
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,656评论 3赞 385
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,322评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,629评论 1赞 293
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,691评论 2赞 312
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,445评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,299评论 3赞 313
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,694评论 3赞 299
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,982评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,244评论 1赞 251
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,642评论 2赞 342
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,829评论 2赞 335