参考来源如何用WGDI进行共线性分析(上) - 简书 (jianshu.com)分析基因组共线性、计算Ks和鉴定WGD —— WGDI | 生信技工 (yanzhongsino.github.io)
还没写完 插个眼
#!/bin/bash
cwd=$(pwd)
bin=$cwd/bin
sample_genome_fa=""
query_genome_fa=""
sample_gff=""
query_gff=""
sample_pep_fa=""
query_pep_fa=""
prefix=""
execute_d_command=false
# 显示脚本的使用方法
show_usage() {
echo "Usage: $0 [-prefix <prefix>] [-s <sample_gff>] [-q <query_gff>] [-sf <sample_genome_fa>] [-qf <query_genome_fa>] [-sp <sample_pep_fa>] [-qp <query_pep_fa>]"
echo "[option]"
echo " -prefix Prefix for the output files"
echo " -s Path to the sample GFF file"
echo " -q Path to the query GFF file"
echo " -sf Path to the sample genome FASTA file"
echo " -qf Path to the query genome FASTA file"
echo " -sp Path to the sample peptide FASTA file"
echo " -qp Path to the query peptide FASTA file"
echo " -d Execute "wgdi -d input.conf" when -d option is provided"
echo " -ks Execute "wgdi -ks input.conf" when -ks option is provided"
echo " --help Show this help message"
}
# 处理命令行选项
while getopts ":prefix:s:q:sf:qf:sp:qp:d-:" opt; do
case $opt in
prefix)
prefix=$OPTARG
;;
s)
sample_gff=$OPTARG
;;
q)
query_gff=$OPTARG
;;
sf)
sample_genome_fa=$OPTARG
;;
qf)
query_genome_fa=$OPTARG
;;
sp)
sample_pep_fa=$OPTARG
;;
qp)
query_pep_fa=$OPTARG
;;
d)
execute_d_command=true
;;
-)
case $OPTARG in
help)
show_usage
exit 0
;;
*)
echo "Invalid option: --$OPTARG" >&2
show_usage
exit 1
;;
esac
;;
\?)
echo "Invalid option: -$OPTARG" >&2
show_usage
exit 1
;;
:)
echo "Option -$OPTARG requires an argument." >&2
show_usage
exit 1
;;
esac
done
# 检查是否提供了prefix参数，如果未提供则使用默认的提取操作
if [ -z "$prefix" ]; then
if [ -z "$query_genome_fa" ]; then
echo "Error: You must provide query_genome.fa parameters." >&2
exit 1
fi
# 提取 sample_genome.fa 和 query_genome.fa 中第一个标点符号前的文字部分作为默认prefix
if [ -z "$sample_gff" ]; then
sample_gff=$query_genome_fa
prefix=$(echo "$sample_genome_fa" | awk -F'.' '{print $1}')_"$(echo "$query_genome_fa" | awk -F'.' '{print $1}')"
fi
if [ -z "$sample_gff" ]; then
sample_gff=$query_genome_fa
mkdir -p $prefix
### make blast database
makeblastdb -in $query_pep_fa -dbtype prot -out $cwd/db/$(echo "$query_pep_fa" | awk -F'.' '{print $1}')
blastp -num_threads 32 -db $cwd/db/$(echo "$query_pep_fa" | awk -F'.' '{print $1}') -query $sample_pep_fa -outfmt 6 -evalue 1e-5 -num_alignments 20 -out $cwd/$prefix/sample.blastp.txt &
python $bin/01.getgff.py $query_gff $prefix/temp.gff
python $bin/02.gff_lens.py $prefix/temp.gff $prefix/input.gff $prefix/input.len
python $bin/03.seq_newname.py $prefix/input.gff $query_pep_fa $prefix/input.pep.fa
python $bin/03.seq_newname.py $prefix/input.gff $query_cds_fa $prefix/input.cds.fa
### or use bin/generate_conf.py
#python $bin/generate_conf.py -p $prefix $sample.genome.fa $sample.gff
else
# 文件路径不同的情况下执行的命令
python $bin/01.getgff.py "$sample_gff" temp_sample.gff
python $bin/01.getgff.py "$query_gff" temp_query.gff
fi
if $execute_d_command; then
gff1= $sample_gff
echo "Running wgdi -d."
echo "[dotplot]
blast = blast file
gff1 = $sample_gff
gff2 = $query_gff
lens1 = lens1 file
lens2 = lens2 file
genome1_name = $(echo "$sample_genome_fa" | awk -F'.' '{print $1}')
genome2_name = $(echo "$query_genome_fa" | awk -F'.' '{print $1}')
multiple = 1
score = 100
evalue = 1e-5
repeat_number = 10
position = order
blast_reverse = false
ancestor_left = ancestor file or none
ancestor_top = ancestor file or none
markersize = 0.5
figsize = 10,10
savefig = savefile(.png, .pdf, .svg)
">input.conf
wgdi -d input.conf
else
echo "Running specific command when sample_gff is not provided."
# 在这里添加您要执行的特定命令（不带 -d 选项）
fi
# 绘制共线性点阵图
wgdi -d ? >input.conf
wgdi -d input.conf
#获取共线性区块
wgdi -icl ? >>input.conf
wgdi -icl input.conf
#计算共线性区块的基因对间的ka和ks
wgdi -ks ? >>input.conf
wgdi -ks ks.total.conf
#整合共线性区块信息
wgdi -bi ? >>input.conf
wgdi -bi input.conf
#过滤并绘制ks点阵图
wgdi -bk ? >>input.conf
wgdi -bk input.conf
#过滤并绘制ks频率分布图
wgdi -kp ? >>input.conf
wgdi -kp input.conf
# 高斯拟合ks频率分布图的峰
wgdi -pf ? >>peak.conf
wgdi -pf peak.conf
#拟合结果作图
wgdi -kf ? >>input.conf
wgdi -kf input.conf
我好菜啊
网友评论