目录
你可以简单地通过conda安装这个软件(AnchorWave还未发布新的版本v1.2.6,因此你可能需要手动从github下载开发版本的AnchorWave,并将二进制文件复制到你的conda环境的bin目录下):
conda create -n quota_Anchor bioconda::quota_anchor
quota_Anchor -h
用法: quota_Anchor [-h] [-v] {longest_pep,longest_cds,pre_col,col,get_chr_length,dotplot,circle,line,ks,class_gene,kde,kf,trios,correct} ...
使用AnchorWave中实现的最长路径算法并考虑基因方向和全基因组复制信息来识别对一对基因组的共线性基因
选项:
-h, --help 显示帮助信息然后退出。
-v, --version 显示版本信息然后退出。
基因共线性分析:
{longest_pep,longest_cds,pre_col,col,get_chr_length,dotplot,circle,line,ks,class_gene,kde,kf,trios,correct}
longest_pep 基于基因组文件和GFF文件调用gffread生成物种蛋白质序列,接着结合GFF文件提取其最长转录本。
longest_cds 基于基因组文件和GFF文件调用gffread产生物种的基因编码序列,接着结合GFF文件提取其最长编码序列。
pre_col 基于染色体长度文件和blast文件产生共线性分析所需的输入文件(称为表文件或包含基因位置信息的blast文件)。
col 基于表文件生成共线性文件。
get_chr_length 根据fai文件和GFF文件产生包含染色体长度和基因总数信息的长度文件。
dotplot 产生共线性基因对点图和同源基因对点图。
circle 共线性结果可视化(圈图)。
line 共线性结果可视化(线形图)。
ks 并行计算共线性基因对的同义和非同义替换率。
class_gene 将基因或者基因对分类为全基因组复制重复、串联重复、近端重复、转座子重复以及离散重复,对于基因分类来说还有不存在同源基因的单基因类别(singleton)。
kde 焦点物种 所有共线性基因对ks/共线性区块中位数ks 的柱形图和高斯核密度评估曲线。
kf 对玉米全基因组复制事件ks进行核密度评估及高斯近似函数拟合并绘图或结合纠正后的物种分化ks峰值绘图。
trios 根据newick格式的二叉树产生trios(焦点物种、姊妹物种以及外群物种组成的三元组)和物种对文件。
correct 拟合物种分化事件的ks峰值并将其矫正到焦点物种的进化速率水平上。
自从玉米与高粱分化以来,玉米相对高粱多经历了一次全基因组复制, 但随后的染色体融合导致这两个物种的染色体数目相同(n = 10)。AnchorWave 最多可以为每个高粱基因提供两条共线性路径,而为每个玉米基因提供一条共线性路径。更多信息请参考文档。
当前工作目录包含 fasta 格式的基因组文件和 gff 格式的基因组注释文件。
wget https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/Zm-B73-REFERENCE-NAM-5.0.fa.gz
wget https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.gff3.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/fasta/sorghum_bicolor/dna/Sorghum_bicolor.Sorghum_bicolor_NCBIv3.dna.toplevel.fa.gz
wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/gff3/sorghum_bicolor/Sorghum_bicolor.Sorghum_bicolor_NCBIv3.57.gff3.gz
gunzip *gz
方便起见,将文件重命令如下。
├── maize.fa
├── maize.gff3
├── sorghum.fa
└── sorghum.gff3
该过程主要包括两个步骤:
- 基于遗传密码表,从基因组文件和注释文件中提取蛋白质序列。
- 对于每个基因,提取其最长转录本。
quota_Anchor longest_pep -f sorghum.fa,maize.fa -g sorghum.gff3,maize.gff3 -p sb.p.fa,zm.p.fa -l sorghum.protein.fa,maize.protein.fa -t 2 --overwrite -merge merged.pep.fa
染色体长度文件除了记录了染色体名字外还记录了染色体长度和基因总数信息,随后可用于共线性分析和绘图。请使用quota_Anchor get_chr_length
查看-s
参数的含义。
quota_Anchor get_chr_length -f sorghum.fa.fai,maize.fa.fai -g sorghum.gff3,maize.gff3 -s 0-9:chr -o sorghum.length.txt,maize.length.txt --overwrite
- 使用 DIAMOND/BLASTP/BLASTN 进行蛋白或者编码序列局部比对。
- 将 BLAST 结果和 GFF 文件信息整合到一个表格文件中。
quota_Anchor pre_col -a diamond -rs sorghum.protein.fa -qs maize.protein.fa -db sorghum.database.diamond -mts 20 -e 1e-10 -b sorghum.maize.diamond -rg sorghum.gff3 -qg maize.gff3 -o sb_zm.table -bs 100 -al 0 -rl sorghum.length.txt -ql maize.length.txt --overwrite
-
通过指定
-r -q
参数, 生成共线性结果。quota_Anchor col -i sb_zm.table -o sb_zm.collinearity -r 2 -q 1 -s 0 --overwrite
-
移除
相对倒位基因对,生成所有
共线性结果(不指定-r -q
)。quota_Anchor col -i sb_zm.table -o sb_zm.collinearity -s 1 -a 1 --overwrite
-
保留
相对倒位基因对, 生成所有
共线性结果(不指定-r -q
)。quota_Anchor col -i sb_zm.table -o sb_zm.collinearity -s 0 -a 1 --overwrite
该过程主要包括两个步骤:
- 从基因组文件和注释文件中提取编码序列。
- 对于每个基因,提取其最长编码序列。
quota_Anchor longest_cds -f sorghum.fa,maize.fa -g sorghum.gff3,maize.gff3 -p sb.cds.fa,zm.cds.fa -l sorghum.cds.fa,maize.cds.fa -t 2 --overwrite -merge merged.cds.fa
quota_Anchor ks -i sb_zm.collinearity -a muscle -p merged.pep.fa -d merged.cds.fa -o sb_zm.ks -t 16 --overwrite
-
同源基因对可视化,使用identity作为图例。
quota_Anchor dotplot -i sb_zm.table -o sb_zm.table.identity.png -r sorghum.length.txt -q maize.length.txt -t order -r_label "Sorghum bicolor" -q_label "Zea mays" -w 1500 -e 1200 -use_identity --overwrite
-
共线性基因对可视化,使用identity作为图例。
quota_Anchor dotplot -i sb_zm.collinearity -o sb_zm.collinearity.identity.png -r sorghum.length.txt -q maize.length.txt -t order -r_label "Sorghum bicolor" -q_label "Zea mays" -w 1500 -e 1200 -use_identity --overwrite
-
共线性基因对可视化,使用ks值作为图例。
quota_Anchor dotplot -i sb_zm.collinearity -o sb_zm.collinearity.ks.png -r sorghum.length.txt -q maize.length.txt -t order -r_label "Sorghum bicolor" -q_label "Zea mays" -w 1500 -e 1200 -ks sb_zm.ks --overwrite
玉米和高粱物种间
quota_Anchor circle -i sb_zm.collinearity -o sb_zm.circle.png -q maize.length.txt -r sorghum.length.txt -rn "Sorghum bicolor" -qn "Zea mays" -cf 9 -sf 9 -rm chr,Chr -fs 14,14 --overwrite
高粱物种内
quota_Anchor circle -i sb_sb.collinearity -o sb_sb.circle.png --overwrite -r ../sorghum.length.txt -q ../sorghum.length.txt -rn "sorghum" -qn "sorghum"
-
两物种共线性基因对可视化
quota_Anchor line -i sb_zm.collinearity -o sb_zm.line.png -l sorghum.length.txt,maize.length.txt -n "Sorghum bicolor,Zea mays" --overwrite
-
三个物种共线性基因对可视化(两次两两比较)
quota_Anchor line -i os_sb.collinearity,sb_sv.collinearity -o os_sb_sv.line.png -l os.length.txt,sb.length.txt,sv.length.txt -n "Oryza sativa, Sorghum bicolor,Zea mays" -rm "chr,Chr" -cf 7 -sf 10 -fs 14,14 --overwrite
该流程参考了DupGen_finder,二者在一些方面有所不同。简单来说,我们在非unique模式下的划分条件更加宽松,而在unique模式下的划分条件更加严格。
-
玉米和玉米的共线性分析
quota_Anchor pre_col -a diamond -rs maize.protein.fa -qs maize.protein.fa -db maize.database.diamond -mts 20 -e 1e-10 -b maize.maize.diamond -rg maize.gff3 -qg maize.gff3 -o zm_zm.table -bs 100 -al 0 -rl maize.length.txt -ql maize.length.txt --overwrite quota_Anchor col -i zm_zm.table -o zm_zm.collinearity -s 0 -m 500 -W 5 -E -0.005 -D 25 -a 1 --overwrite
-
下载Musa balbisiana基因组文件和注释文件并重命名
wget https://api.ncbi.nlm.nih.gov/datasets/v2/genome/accession/GCA_004837865.1/download\?include_annotation_type\=GENOME_FASTA\&include_annotation_type\=GENOME_GFF
├── banana.B.fa └── banana.B.gff
-
野蕉和玉米的共线性分析
quota_Anchor longest_pep -f banana.B.fa -g banana.B.gff -p B.p.pep -l banana.B.pep -t 1 --overwrite quota_Anchor get_chr_length -f banana.B.fa.fai -g banana.B.gff -s CM01 -o banana.B.length.txt --overwrite quota_Anchor pre_col -a diamond -rs banana.B.pep -qs maize.protein.fa -db banana.B.database.diamond -mts 20 -e 1e-10 -b banana.B.maize.diamond -rg banana.B.gff -qg maize.gff3 -o bananaB_zm.table -bs 100 -al 0 -rl banana.B.length.txt -ql maize.length.txt --overwrite quota_Anchor col -i bananaB_zm.table -o bananaB_zm.collinearity -s 0 --overwrite -D 25 -a 1
-
玉米基因及同源基因对的划分 Unique模式
quota_Anchor class_gene -b maize.maize.diamond -g maize.gff3 -q zm_zm.collinearity -qr bananaB_zm.collinearity -o maize_classify_dir -p maize -s 1 -d 10 --overwrite -u
Non-unique模式
quota_Anchor class_gene -b maize.maize.diamond -g maize.gff3 -q zm_zm.collinearity -qr bananaB_zm.collinearity -o maize_classify_dir -p maize -s 1 -d 10 --overwrite
这个流程参考了 ksrates, 两者在一些地方有所不同。简单来说,该流程使用基于-r_value -q_value
参数获得的共线性基因对ks值拟合结果作为物种分化峰,而ksrates使用RBH基因对ks值拟合结果作为物种分化峰。此外,拟合方法也有所不同。
以下是当前目录信息以及水稻和狗尾草的数据来源。
wget https://ftp.ebi.ac.uk/ensemblgenomes/pub/release-59/plants/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna_rm.toplevel.fa.gz
wget https://ftp.ebi.ac.uk/ensemblgenomes/pub/release-59/plants/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.59.gff3.gz
wget https://ftp.ebi.ac.uk/ensemblgenomes/pub/release-59/plants/fasta/setaria_viridis/dna/Setaria_viridis.Setaria_viridis_v2.0.dna.toplevel.fa.gz
wget https://ftp.ebi.ac.uk/ensemblgenomes/pub/release-59/plants/gff3/setaria_viridis/Setaria_viridis.Setaria_viridis_v2.0.59.gff3.gz
├── raw_data
│ ├── maize.fa
│ ├── maize.gff3
│ ├── oryza.fa
│ ├── oryza.gff3
│ ├── setaria.fa
│ ├── setaria.gff3
│ ├── sorghum.fa
│ └── sorghum.gff3
└── scripts
├── ks_pipeline.py
└── longest_pipeline.py
-
对于
raw_data(input_dir)
的每一个物种生成其最长转录本和最长编码序列。python ./scripts/longest_pipeline.py -i raw_data -o output_dir --overwrite
-
获得每个物种的染色体长度文件。 你可能需要根据
quota_Anchor get_chr_length
来了解-s
参数的含义。a)
find ./raw_data/*fai |awk '{printf "%s,", $1}' find ./raw_data/*gff3 |awk '{printf "%s,", $1}' find ./raw_data/*gff3 |awk '{printf "%s,", $1}'|sed s/gff3/length\.txt/g
quota_Anchor get_chr_length -f ./raw_data/maize.fa.fai,./raw_data/oryza.fa.fai,./raw_data/setaria.fa.fai,./raw_data/sorghum.fa.fai -g ./raw_data/maize.gff3,./raw_data/oryza.gff3,./raw_data/setaria.gff3,./raw_data/sorghum.gff3 -s 0-9,CHR,chr,Chr:0-9,CHR,chr,Chr:0-9,CHR,chr,Chr:0-9,CHR,chr,Chr -o ./raw_data/maize.length.txt,./raw_data/oryza.length.txt,./raw_data/setaria.length.txt,./raw_data/sorghum.length.txt --overwrite
b)
quota_Anchor get_chr_length -f "$(find ./raw_data/*fai |awk '{printf "%s,", $1}')" -g "$(find ./raw_data/*gff3 |awk '{printf "%s,", $1}')" -s 0-9,CHR,chr,Chr:0-9,CHR,chr,Chr:0-9,CHR,chr,Chr:0-9,CHR,chr,Chr -o "$(find ./raw_data/*gff3 |awk '{printf "%s,", $1}'|sed s/gff3/length\.txt/g)" --overwrite
-
根据提供的newick格式的二叉树获得物种对文件和trios文件。
quota_Anchor trios -n "(((maize, sorghum), setaria), oryza);" -k "maize" -ot ortholog_trios_maize.csv -op species_pairs.csv -t tree.txt --overwrite
Species_1 Species_2 q_value r_value get_all_collinear_pairs maize setaria 1 1 0 maize setaria 1 1 0 sorghum setaria 1 1 0 maize oryza 1 1 0 sorghum oryza 1 1 0 setaria oryza 1 1 0 -
对于每一个物种对产生其共线性结果及共线性基因对所对应的ks值。 注:
./scripts/ks_pipeline.py
这个脚本在共线性过程使用species_pairs.csv
的Species_1
列的值作为查询物种,使用Species_2`列的值作为参考物种。./scripts/ks_pipeline.py
脚本会根据species_pairs.csv
物种对文件的q_value
,r_value
和get_all_collinear_pairs
调整共线性过程的参数。- 你可能需要根据
quota_Anchor col
来了解这三个参数的含义或者参考文档。
python ./scripts/ks_pipeline.py -i raw_data -o output_dir -s species_pairs.csv -a diamond -l raw_data --overwrite -plot_table
-
对于每个分化峰进行拟合,并且基于trios三元组将速率矫正到所关注物种玉米的进化速率上。 注:
find ./output_dir/02synteny/*0.ks |awk '{printf "%s,", $1}'
命令中的0
表示物种对文件的get_all_collinear_pairs
列的值。- 物种对文件(由-s参数指定,species_pairs.csv)中物种对的顺序必须与ks文件(由-k参数指定)的顺序一致
find ./output_dir/02synteny/*0.ks |awk '{printf "%s,", $1}'
quota_Anchor correct -k "./output_dir/02synteny/maize_sorghum0.ks,./output_dir/02synteny/maize_setaria0.ks,./output_dir/02synteny/sorghum_setaria0.ks,./output_dir/02synteny/maize_oryza0.ks,./output_dir/02synteny/sorghum_oryza0.ks,./output_dir/02synteny/setaria_oryza0.ks" -s species_pairs.csv -t ortholog_trios_maize.csv -kr 0,2 -ot outfile_divergent_peaks.csv --overwrite
-
玉米全基因组复制事件共线性基因对ks值柱形图及高斯核密度评估曲线。
quota_Anchor pre_col -a diamond -rs ./output_dir/01longest/maize.longest.pep -qs ./output_dir/01longest/maize.longest.pep -db ./maize/maize.database.diamond -mts 20 -e 1e-10 -b ./maize/maize.maize.diamond -rg ./raw_data/maize.gff3 -qg ./raw_data/maize.gff3 -o ./maize/zm_zm.table -bs 100 -al 0 --overwrite quota_Anchor dotplot -i ./maize/zm_zm.table -o ./maize/zm.zm.png -r ./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -r_label maize -q_label maize -use_identity --overwrite quota_Anchor col -i ./maize/zm_zm.table -o ./maize/zm_zm.collinearity -a 1 -m 1000 -W 1 -D 25 -I 5 -E -0.01 -f 0 --overwrite quota_Anchor dotplot -i ./maize/zm_zm.collinearity -o ./maize/zm.zm.collinearity.png -r ./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -r_label maize -q_label maize -use_identity --overwrite quota_Anchor ks -i ./maize/zm_zm.collinearity -a mafft -p ./output_dir/01longest/maize.longest.pep -d ./output_dir/01longest/maize.longest.cds -o ./maize/zm.zm.ks -t 16 --overwrite quota_Anchor dotplot -i ./maize/zm_zm.collinearity -o ./maize/zm.zm.collinearity.ks.png -r ./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -r_label maize -q_label maize --overwrite -ks ./maize/zm.zm.ks
quota_Anchor kde -i ./maize/zm_zm.collinearity -r./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -o ./maize/zm.zm.kde.png -k ./maize/zm.zm.ks --overwrite
你需要根据
zm.zm.table.png
或者其他方法提供玉米wgd峰值的数量。quota_Anchor kf -i ./maize/zm_zm.collinearity -r ./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -o ./maize/zm.zm.kf.png --overwrite -k ./maize/zm.zm.ks -components 2 -f maize -kr 0,2
-
利用高斯混合模型对wgd基因对ks进行分组,对各组分进行核密度估计和高斯近似拟合,取382次核密度估计的众数为最初的物种分化峰,并基于三重奏将分化峰校正至焦点物种水平。
quota_Anchor kf -i ./maize/zm_zm.collinearity -r ./raw_data/maize.length.txt -q ./raw_data/maize.length.txt -o ./maize/zm.zm.png --overwrite -k ./maize/zm.zm.ks -components 2 -f maize -correct_file outfile_divergent_peaks.csv -kr 0,2
常见问题请看FAQ_zh.md。