基因组直系同源基因kaks和LTR插入时间计算

基因组直系同源基因kaks计算

#只保留第一个和第二个空格之间的内容作为序列名，用于下载的多个基因组的数据格式化

seqkit replace -p '^.+?\s([^\s]+)\s.*' -r '$1' input.fasta > output.fasta
#修改所有蛋白序列名称，只留下 "> 基因名"，删掉后面的注释信息。
sed -E 's/^>([^[:space:]]+).*/>\1/' CY.pep > CY1.pep
# 把两个物种的蛋白序列放在一个文件夹下
orthofinder -f LmigXrip/
# 运行后会显示职别的同源基因数目、蛋白分配比例、正交群数目和单拷贝基因数目。
# 合并单拷贝基因的蛋白序列文件即可
cat *fa > 2.pep
#整理蛋白和cds序列名，去除空格之后的其他内容，防止axt格式出错
seqkit replace -p " .*" -r "" input.fasta > output.fasta
#取出基因对文件，针对两个物种分别取出子文件
seqkit fx2tab all.pep > 2pep

###############脚本按照2pep的蛋白名提取cds####################
#!/bin/bash

# 用法: ./prepare_paraat_input.sh <cds_file> <output_dir>
set -e

CDS_FILE="$1"
OUT_DIR="$2"

# 检查输入参数
if [ ! -f "$CDS_FILE" ]; then
  echo "错误：CDS 文件不存在！"
  exit 1
fi

if [ ! -d "$OUT_DIR" ]; then
  echo "创建输出目录：$OUT_DIR"
  mkdir -p "$OUT_DIR"
fi

cd "$OUT_DIR"

echo "Step 1: 检查 all.pep 文件是否存在"
if [ ! -f "all.pep" ]; then
  echo "错误：all.pep 文件不存在，请确保 all.pep 文件位于输出目录中！"
  exit 1
fi

echo "Step 2: 整理蛋白和CDS序列名（去除空格后内容）"
seqkit replace -p " .*" -r "" all.pep > 2.pep
seqkit replace -p " .*" -r "" "$CDS_FILE" > 2.cds

echo "Step 3: 提取蛋白ID并拆分为两组"
seqkit fx2tab 2.pep > 2pep.tab
cut -f 1 2pep.tab > 2_id

awk 'NR%2==1' 2_id > ID1
awk 'NR%2==0' 2_id > ID2
paste ID1 ID2 > 2.id

echo "Step 4: 提取对应CDS并清理名称"
seqkit grep -f ID1 2.cds > cds1.fa
seqkit grep -f ID2 2.cds > cds2.fa
cat cds1.fa cds2.fa > merged.cds
seqkit replace -p " .*" -r "" merged.cds > 2.cds1

echo "准备完成 ✅"
echo "生成文件: 2.pep  2.cds1  2.id"

rm 2.cds
rm merged.cds


#####################

#现在有2.pep,2.cds,2.id ，三个文件可作为ParaAT的输入文件
ParaAT2.pl -h 2.id -n 2.cds1 -a 2.pep -p proc -m mafft -f axt -g -k -o result_dir 
####################################
-h 基因对文件
-n 基因对的核酸序列
-a 基因对的蛋白序列
-m 比对工具，可以用muscle或mafft
-p proc文件必须与输出位置在同一个目录下，这一点非常重要
export PATH=$PATH:/pwd
chmod +x Epal2nal.pl
####################################

LTR插入时间计算

使用EDTA可以获得从头预测的LTR的插入时间，不过使用的是水稻单子叶植物的碱基突变速率，需要换算成自己对应类群的碱基替换速率。我做的苔藓，使用小立碗藓的速率计算。

singularity exec --bind /mnt/d/bio_data/001.repeat:/mnt/d/bio_data/001.repeat /mnt/d/bio_soft/EDTA.sif EDTA_raw.pl \
--genome /mnt/d/bio_data/001.repeat/cy.new.fasta \
--type ltr \
--species others \
--threads 14 \
--overwrite 1 \
--u 1.3e-8
如下路径文件最后一列为插入时间，是使用水稻的碱基突变速率计算得来的
(base) xfm@XFM-R9KP:/mnt/d/bio_data/001.repeat$ ll cy.new.fasta.mod.EDTA.raw/LTR/cy.new.fasta.mod.pass.list
-rwxrwxrwx 1 xfm xfm 19556 Jun 25 15:51 cy.new.fasta.mod.EDTA.raw/LTR/cy.new.fasta.mod.pass.list*
####使用以下脚本修改结果，使用小立碗藓的碱基突变速率，9.4e-09
#!/usr/bin/env python

import sys

if len(sys.argv) != 3:
    print("Usage: {} input.txt mutation_rate".format(sys.argv[0]))
    sys.exit(1)

input_file = sys.argv[1]
mutation_rate = float(sys.argv[2])
output_file = input_file + "2"

default_rate = 1.3e-8

rate_ratio = mutation_rate / default_rate

with open(input_file, "r") as fin, open(output_file, "w") as fout:
    for line in fin:
        if not line.strip() or line.startswith("#"):
            fout.write(line)
            continue

        items = line.strip().split()

        old_time = float(items[-1])
        new_time = old_time / rate_ratio

        items[-1] = str(new_time)

        fout.write("\t".join(items) + "\n")
        
        
#########################