python 自建基因比对数据库基因组比对

转载

云端小梦 2024-01-30 10:00:36

文章标签 python 自建基因比对数据库 linux 服务器动态规划 git 文章分类 Python 后端开发

本文采用两个不同品种的拟南芥进行全基因组比对和变异检测。这种比对方法使得每个相对应的染色体名称都一样。并且对于两个相同物种之间存在倒位等染色体变异，它的全基因组比对过程也是类似的。两个基因组进行从头到尾的碱基水平上的全基因组比对。

1.下载基因组的序列文件和参考基因组的注释文件

使用gean工具进行从sdi格式转化为fasta格式。

#如果没有安装GEAN，可以通过以下方式进行安装，这只是文件格式的转化，
如果直接下载的格式是fasta格式这个过程是不必要的。
git clone https://github.com/baoxingsong/GEAN.git
cd GEAN
cmake CMakeLists.txt
make

wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr1.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr2.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr3.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr4.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr5.fas
#上文是下载的每条染色体，使用cat命令进行合成一个。
cat TAIR10_chr1.fas TAIR10_chr2.fas TAIR10_chr3.fas TAIR10_chr4.fas TAIR10_chr5.fas > tair10.fa

wget http://mtweb.cs.ucl.ac.uk/mus/www/19genomes/variants.SDI/bur_0.v7c.sdi
gean pseudogeno -r tair10.fa -v bur_0.v7c.sdi -o bur_0.fa

如果得到的两个基因组的名称不一样则需要进行，名称规范化。

2. 提取参考基因组的中的保守序列和将保守序列分别匹配到两个基因组上

-x splice 代表将保守序列映射到基因组，-t 代表线程数，-k 代表kmer， -a代表输出sam格式的文件，-p和-N代表如果它们的比对得分大于初次比对得分的0.4将会保留最多保留20个结果，还要输入他们的基因组和保守序列。

anchorwave gff2seq -r tair10.fa -i TAIR10_GFF3_genes.gff -o cds.fa

minimap2 -x splice -t 6 -k 12 -a -p 0.4 -N 20 tair10.fa cds.fa > ref.sam
minimap2 -x splice -t 6 -k 12 -a -p 0.4 -N 20 bur_0.fa cds.fa > bur_0.sam

3.对两个品种拟南芥基因组进行变异检测和全基因组比对

在这个例子中我们使用“genoAli”命令而不是上一篇文章中的“proali”-v输出变异检测结果的vcf文件，-i代表输入参考基因组的注释文件，-r 和-s 参考和查询基因组文件，-o和-f 输出变异比对结果的maf格式文件，-n文件可以用于可视化结果。如果过程中消耗的内存过大，可以减小滑动窗口的大小。此外，-m设置最小的外显子长度，不要随意修改，并且和上个步骤一样。

anchorwave genoAli -i TAIR10_GFF3_genes.gff -as cds.fa -r tair10.fa -a bur_0.sam -ar ref.sam -s bur_0.fa -v bur_0.vcf -n bur_0.anchors -o bur_0.maf -f bur_0.f.maf -t 1

对输出的共线性进行作图，代码如下

library(ggplot2)
#设置一个坐标函数
changetoM <- function ( position ){
  position=position/1000000;
  paste(position, "M", sep="")
}
#读取数据
data =read.table("bur_0.anchors", header=TRUE)
#筛选染色体并设置为因子
data = data[which(data$refChr %in% c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5")),]
data = data[which(data$queryChr %in% c("Chr1", "Chr2", "Chr3","Chr4","Chr5")),]
data$refChr = factor(data$refChr, levels=c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5"))
data$queryCh = factor(data$queryChr, levels=c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5" ))
#使用ggplot2包画图并美化它.
ggplot(data=data, aes(x=queryStart, y=referenceStart))+
  geom_point(size=0.5, aes(color=strand)) +
  facet_grid(refChr~queryChr, scales="free", space="free") +
  labs(x="Query_bur_0", y="Reference_tair10")+scale_x_continuous(labels=changetoM) +
  scale_y_continuous(labels=changetoM) +
  theme(axis.line = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(fill =NA,color="black", size=0.5, linetype="solid"),
        axis.text.y = element_text( colour = "black"),
        legend.position='none',
        axis.text.x = element_text(angle=300,hjust=0, vjust=0.5, colour = "black") )

python 自建基因比对数据库基因组比对_python 自建基因比对数据库