JCVI 基因组共线性分析与可视化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

GENOME_PATH=/path/to/genome/fasta

function formated() {
    name=$(basename ${1%%.*})  ## 输入 gff3 文件，提取物种名
    mkdir -p temp_dir

    # GFF3 转 BED 格式
    # --primary_only 只保留每个基因的一个转录本（去除冗余 isoform）
    python -m jcvi.formats.gff bed --type=mRNA --key=ID --primary_only $1 -o ${name}.bed

    # 提取 CDS 序列
    gffread -g ${GENOME_PATH}/${name}.fa -x ./temp_dir/${name}.cds $1

    # 去重：基于 BED 文件去除重复基因
    python -m jcvi.formats.bed uniq ${name}.bed
    mv ${name}* ./temp_dir

    # 根据去重后的 BED 提取对应的 CDS 序列
    seqkit grep -f <(cut -f 4 ./temp_dir/${name}.uniq.bed) ./temp_dir/${name}.cds \
      | seqkit seq -i > ./temp_dir/${name}.uniq.cds

    # 创建符号链接（JCVI 默认在当前目录查找）
    ln -s ./temp_dir/${name}.uniq.cds ./${name}.cds
    ln -s ./temp_dir/${name}.uniq.bed ./${name}.bed
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

function pairwiseSynteny() {
    query=${1}
    target=${2}
    if [ "${query}" != "${target}" ]; then
        # 计算正交同源基因（Ortholog）
        python -m jcvi.compara.catalog ortholog ${query} ${target} --full --no_strip_names

        # 过滤 anchors，保留最小跨度 ≥ 30 的共线性区块
        python -m jcvi.compara.synteny screen --minspan=30 --simple \
            ${query}.${target}.anchors ${query}.${target}.anchors.new
    fi
}

# 批量运行：从 compares 文件中读取物种对
cat ./compares | while read a b
do
    pairwiseSynteny $a $b
    echo "Finish $a & $b!"
done

1
2
3
4
5

# 合并所有物种的 BED 文件
python -m jcvi.formats.bed merge \
    Oryza_sativa.bed Oryza_rufipogon.bed Oryza_indica.bed \
    Oryza_barthii.bed Oryza_nivara.bed Oryza_glaberrima.bed \
    -o all.bed

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 将各 pairwise 的 RBH（Reciprocal Best Hit）结果合并
python -m jcvi.formats.base join \
    Oryza_rufipogon.Oryza_sativa.rbh.simple \
    Oryza_indica.Oryza_sativa.rbh.simple \
    Oryza_barthii.Oryza_sativa.rbh.simple \
    Oryza_nivara.Oryza_sativa.rbh.simple \
    Oryza_glaberrima.Oryza_sativa.rbh.simple \
    --column=1,1,1,1,1 \
    -o merged.rbh

# 重排列顺序（将参考物种放在第一列）
cat merged.rbh | awk -v OFS="\t" '{print $2,$1,$3,$5,$7,$9}' > all.blocks

1
2
3

python -m jcvi.graphics.synteny \
    all.blocks all.bed all.layout \
    --glyphcolor=orthogroup --glyphstyle=arrow