搜索
查看: 12052|回复: 2

使用chromPlot画染色体图

[复制链接]

3

主题

17

帖子

380

积分

中级会员

Rank: 3Rank: 3

积分
380
发表于 2017-6-24 13:32:15 | 显示全部楼层 |阅读模式
circos的图都有视觉审美疲劳了,这里推荐一个R包chromPlot。
如何使用请自己看manual。以下是我用棉花的基因组数据使用chromPlot画棉花染色体基因,差异基因和转座子的分布来介绍下这个R包

效果图 Fig1

1 使用chromPlot在画布上画出棉花的26条染色体
1.1 获得每条染色体的长度,使用大神开发的seqkit,原谅我方法比较糙
[AppleScript] 纯文本查看 复制代码
wget [url]https://github.com/shenwei356/seqkit/releases[/url]
tar -zxvf *.tar.gz
seqkit seq -wo cotton.fa > Ghgenome.fasta # 
for i in {1..9}; do  grep -A1 ">A0$i" Ghgenome.fasta > A${i}.fa;done;
for i in {1..9}; do  grep -A1 ">D0$i" Ghgenome.fasta > D${i}.fa;done;
for i in {10..13}; do  grep -A1 ">A$i" Ghgenome.fasta > A${i}.fa;done;
for i in {10..13}; do  grep -A1 ">D$i" Ghgenome.fasta > A${i}.fa;done;
../../../seqkit stats * |awk '{print $1"\t"1"\t"$5"\t"$1}'|perl -pi -e 's/\.fa//g'


1.2获得染色体,并输出我们想要的格式
棉花是异源多倍体,有2个基因组A和D,分别以A01-A13, D01-D13表示
[AppleScript] 纯文本查看 复制代码
 ../../../seqkit stats * |awk '{print $1"\t"1"\t"$5"\t"$1}'|perl -pi -e 's/\.fa//g'|perl -pi -e 's/A/A0/g if /A\d{1}\t/ '| perl -pi -e 's/D/D0/g if /D\d{1}\t/'


整理后的数据格式如下
[AppleScript] 纯文本查看 复制代码
Chrom	Start	End	Name
A10	1	100,866,604	A10
A11	1	93,316,192	A11
A12	1	87,484,866	A12
A13	1	79,961,121	A13
A01	1	99,884,700	A01

1.3 R代码画染色体
[AppleScript] 纯文本查看 复制代码
library(chromPlot)
setwd("/Users/zt/Desktop/chrplot")
table=read.table("chrom.txt",header = T,sep = "\t",stringsAsFactors = F)
table
chromPlot(gaps=table)


Fig2 画染色体

2 对图进行升级,在染色体上画出转座子信息,棉花的转座子注释信息我提前储存在Ghte_sort.bed
棉花中Gypsy在染色体的转座子比例中占了绝大部分,我用蓝色去表示它。其他转座子用红色,仍旧整理成5列的数据,通过这张图我们可以看出,Gypsy在染色体上的分布
2.1 整理要读取的转座子注释文件,清除位于scaffold的垃圾数据
[AppleScript] 纯文本查看 复制代码
[/b]cut -f 1,2,3,4 Ghte_sort.bed|perl -ne 'print unless /scaff|Novel/' |awk '{print $0"\tblue"}' |perl -pi -e 's/blue/red/g unless /Gypsy/ ' > gene_annotation.bed


整理后的数据格式如下
[AppleScript] 纯文本查看 复制代码
Chrom	Start	End	Name	Colors
A01	5	334	LTR/Gypsy	blue
A01	14191	14632	Unknown	red
A01	21929	22035	MobileElement	red
A01	23421	23605	Unknown	red
A01	23492	23635	Unknown	red


2.3 画包含转座子信息的染色体, Fig. 3
[AppleScript] 纯文本查看 复制代码
te_annotation=read.table("te_annotation.bed",header = T,sep = "\t",stringsAsFactors = F)
chromPlot(gaps=table, bands=te_annotation,annot1 = gene, segment = DEgene,figCols=3,chr = c("A01","D01"))


3 添加基因和差异基因的信息3.1数据准备
circos中我们经常可以看到一圈圈的基因密度,CG含量,如何 把基因的密度画上去,比如左边画基因的密度,右边画差异基因,可以很好的体现分布的规律#基因的格式文件如下
[AppleScript] 纯文本查看 复制代码
Chrom	Start	End	Name
A01	15705	19194	Gh_A01G0001
A01	22808	24529	Gh_A01G0002
A01	36428	36860	Gh_A01G0003
A01	40962	41405	Gh_A01G0004
A01	48194	49051	Gh_A01G0005
A01	55552	59030	Gh_A01G0006

[AppleScript] 纯文本查看 复制代码
# 差异基因的格式如下
Chrom	Start	End	Group
A01	9369	14076	DGE
A01	161390	161756	DGE
A01	814253	815673	DGE
A01	1125455	1126876	DGE


3.2 把差异基因的两列以柱形图形式画上去,以下为Fig. 4 R代码
[AppleScript] 纯文本查看 复制代码
gene= read.table("gene_2.txt",header = T,sep = "\t")
DEgene=read.table("DGE2_gene.txt",header=T,sep="\t")
DEgene                
chromPlot(gaps=table, bands=gene_annotation,annot1 = gene, segment = DEgene,figCols=3,chr = c("A01","D01"))


4 添加分子标记4.1 分子标记信息整理
最常见的工作是把分子标记或者基因家族在染色体标上去分子标记个的数据格式如下
[AppleScript] 纯文本查看 复制代码
Chrom     Start       End       ID    Colors
A01	95813373	95814820	XLOC_013344	black
A01	97199488	97199690	XLOC_028222	black

4.2 以下是R代码4.2 以下是R代码
[AppleScript] 纯文本查看 复制代码
 
mark=read.table("marker.txt",header = T)
chromPlot(gaps=table,stat = mark, statCol="Value", statName="Value",bands=gene_annotation,annot1 = gene, segment = DEgene,figCols=3,chr = c("A01","D01"))
​









以上是我学习该包的流程,向大家介绍这个包






Fig1

Fig1

Fig2

Fig2

Fig3

Fig3

Fig 4

Fig 4
回复

使用道具 举报

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-6-25 19:41:33 | 显示全部楼层
获得每条染色体的长度的方法的确比较糙
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

1

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2022-5-19 16:44:57 | 显示全部楼层
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2023-1-29 23:46 , Processed in 0.157027 second(s), 38 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.