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