litover只能提供一些热门基因组的转换，而西农开发的RGD只有一小部分基因组，不过非常棒的是他们提供了maf文件，这是手动创建chain文件中最耗时间和资源的一步。本文主要记录如何以西农提供的文件获取chain文件并最终做到坐标的转换。
当然如果RGD没有提供所需的chain文件，也可以参照他们的脚本自主创建，不过这一步需要耗费大量的时间和服务器资源（sheep1比对cattle1花费了我 24个小时和近200g的内存）

文件下载和准备

1 maf文件

http://animal.omics.pro/code/index.php/RGD
非常棒的数据库，有很多有用的工具和文件

maf文件下载：RGD显示 (omics.pro)
下载完毕后建议用md5校验一下：

md5sum -c md5.txt 

该论文的附表中提供了上述maf文件的参考基因组信息：RGD v2.0：反刍动物功能和进化基因组学数据库的重大更新 |核酸研究 |牛津学术 (oup.com)

脚本位置：RGD显示 (omics.pro)

2 软件下载

last： Martin Frith / last · GitLab（下载很奇怪）
gtfToGenePred：https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/gtfToGenePred
faToTwoBit：https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
axtChain：wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtChain
chainMergeSort：wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainMergeSort
liftover：https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver
multiz：GitHub - multiz/multiz：DNA 多序列比对器，宾夕法尼亚州立大学米勒实验室的正式版本（直接conda）（不对，好像用不上，不改了）
这几个ucsc的软件下载后修改一下权限可以直接用

其余使用但没有提到的程序都包含在last内。

3 参考基因组文件

1.RGD提供的maf文件中包含"NW"染色体，需要从NCBI下载参考科基因组（可以使用RGD提供的下载链接RGD显示 (omics.pro)
）
2.常染色体和性染色体是数字形式，需要进行染色体号的转换

awk 'NR==FNR{a[$1]=$2;next}{for(i in a){sub(i,a[i])};print}' replace.txt input.file > out.file
代码取自：https://www.cnblogs.com/irockcode/p/7906192.html

后续又发现一个问题，maf文件中的染色体号并不是标准的数字染色体号，而是 Sheep.1 ARS_UCD_addY.1 这种形式，
为了便于操作，直接使用了下述代码
sed 's/>/>ARS_UCD_addY./g'  GCF_002263795.1_ARS-UCD1.2_genomic.chr.addY.fna  > ARS_UCD_addY.GCF_002263795.1_ARS-UCD1.2_genomic.chr.addY.fna 

3.RGD的参考基因组均进行了加Y处理。其中牛的是5.1，羊的不清楚（推测可能是ASM1117029v1这个由西农自主组装的基因组）
ASM1117029v1：绵羊基因组组装 ASM1117029v1 - NCBI - NLM (nih.gov)
Btau_5.0.1：Bos taurus 基因组组装 Btau_5.0.1 - NCBI - NLM (nih.gov)
别忘了下载对应的注释文件

参考基因组的转换和合并

pip install pyfaidx #别用conda下载
faidx -x genome.fa #拆分参考基因组，只要其中的Y染色体
sed -i 's/CM022046.1/Y/g' CM022046.1.fna
sed -i 's/NC_016145.1/Y/g' NC_016145.1.fna
grep '>' GCF_002263795.1_ARS-UCD1.2_genomic.fna > tt
awk '/NW_020190115.1/ {print NR}' GCF_002263795.1_ARS-UCD1.2_genomic.fna #这两步是为了确定第一个非常染色体和性染色体的行号，如NW_020190115.1是32854981
sed '32854980 r NC_016145.1.fna' GCF_002263795.1_ARS-UCD1.2_genomic.chr.fna > output_file #希望没错，主要是这也没法检查，还好我用不到X,Y

修改注释文件
牛只有gff文件

conda install gffread #我更喜欢用gtf，所以把提供的gff转为gtf
gffread GCF_000003205.7_Btau_5.0.1_genomic.gff -T -o GCF_000003205.7_Btau_5.0.1_genomic.gtf
grep "NC_016145.1"  GCF_000003205.7_Btau_5.0.1_genomic.gtf > NC_016145.1.gtf
sed -i 's/NC_016145.1/Y/g' NC_016145.1.gtf
awk 'NR==FNR{a[$1]=$2;next}{for(i in a){sub(i,a[i])};print}' NC-chr-cattle.txt GCF_002263795.1_ARS-UCD1.2.gtf > GCF_002263795.1_ARS-UCD1.2.chr.gtf 
cat  GCF_002263795.1_ARS-UCD1.2.chr.gtf  NC_016145.1.gtf > GCF_002263795.1_ARS-UCD1.2.chr.addY.gtf
同理，增加一行改变染色体号的命令
grep -o '^[^#]*'  GCF_002263795.1_ARS-UCD1.2.chr.addY.gtf > GCF_002263795.1_ARS-UCD1.2.chr.addY.#.gtf  #先删除注释行
sed 's/^/ARS_UCD_addY./' GCF_002263795.1_ARS-UCD1.2.chr.addY.#.gtf > ARS_UCD_addY.GCF_002263795.1_ARS-UCD1.2.chr.addY.gtf

绵羊只有gbff注释文件

找到的一些转换脚本都不是很好用，先这样把，不要y染色体注释文件了

流程

分离maf文件

从RGD下载的maf文件是以某一个参考基因组为基准，包含110个物种的比对信息，因此首先要做的是从文件中分离出需要的物种。
谢谢你，Newbing

#!/bin/bash

# 检查参数是否正确
if [ $# -ne 2 ]; then
    echo "用法: $0 <输入文件> <输出文件>"
    exit 1
fi

input_file="$1"
output_file="$2"

# 读取输入文件的每一行
while IFS= read -r line; do
    # 如果是注释行（以#开头），直接输出到输出文件
    if [[ "$line" == "#"* ]]; then
        echo "$line" >> "$output_file"
    elif [[ "$line" == "a"* ]]; then
        # 如果是以a开头的行，开始一个新区块
        block=""
        block="$block$line"$'\n'
        # 读取下一行，直到遇到下一个以a开头的行或文件结束
        while IFS= read -r next_line; do
             if [ -z "$next_line" ]; then
            break
            fi
            block="$block$next_line"$'\n'
        done
        # 判断区块中是否有以s Sheep开头的行
        if [[ "$block" == *"s Sheep"* ]]; then
            # 输出区块的前两行和以s Sheep开头的行
            echo "$block" | head -n 2 >> "$output_file"
            echo "$block" | grep "s Sheep" >> "$output_file" #只需修改这里的grep内容即可筛选不同的参考基因组
            echo "" >> "$output_file"  # 输出一行空行
        fi
    fi
done < "$input_file"

bash tiqu.sh chr1.maf chr1.sheep.maf#这里整个循环挺不错的

for i in chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chr23 chr24 chr25 chr26 chr27 chr28 chr29 chrX chrY
do
gzip -d ${i}.maf.gz
bash tiqu.sh ${i}.maf ${i}.sheep.maf
head -1 chrY.maf > head #为下一步做准备，主要是为了提取个注释行，预计会有30个报错（head: 无法打开'chrY.maf' 读取数据: 没有那个文件或目录）
gzip "${i}".maf
done
全部解压缩的话会非常占硬盘空间（有钱人请无视重新压缩那步）

运行起来很慢，建议手动多线程

合并maf文件

这里应该能有更完美的处理方法，可惜我不会
for i in "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 26 27 28 29 X Y"
do
cat head chr${i}.sheep.maf >> head
done

另一方面，ucsc也有合并chain文件的软件，见后文

psl文件

maf-convert psl all.sheep.maf > all.sheep.maf.psl

处理参考基因组和注释文件

./faToTwoBit GCF_002742125.1_Oar_rambouillet_v1.0_genomic.chr.addY.fa GCF_002742125.1_Oar_rambouillet_v1.0_genomic.chr.addY.fa.2bit
./faToTwoBit GCF_002263795.1_ARS-UCD1.2_genomic.chr.addY.fna GCF_002263795.1_ARS-UCD1.2_genomic.chr.addY.fna.2bit
./gtfToGenePred  -genePredExt -ignoreGroupsWithoutExons -allErrors ARS_UCD_addY.GCF_002263795.1_ARS-UCD1.2.chr.addY.gtf ARS_UCD_addY.GCF_002263795.1_ARS-UCD1.2.chr.addY.gtf.genePred 
注意！！ 对于NCBI下载的注释文件， -allErrors参数是必须的
gtfToGenePre会报一个非常奇怪的错误 
invalid gffGroup detected on line: 10   Gnomon  start_codon     22347475        22347477        0.000000        +       0       gene_id "LOC101908096"; transcript_id "unknown_transcript_1"; 
GFF/GTF group unknown_transcript_1 on 4+, this line is on 10+, all group members must be on same seq and strand
合理怀疑是把unknown_transcript_1当成一个组了，暂时不知道怎么解决这个问题，只能先忽略

chain

./axtChain -linearGap=loose -psl chr1.sheep.maf.middle.psl ARS_UCD_addY.GCF_002263795.1_ARS-UCD1.2_genomic.chr.addY.fna.2bit sheep.GCF_002742125.1_Oar_rambouillet_v1.0_genomic.chr.addY.2bit ARS-UCD1.2ToOar_rambouillet_v1.0.chain #不懂参数，直接复制来的

chain合并

chainMergeSort - Combine sorted files into larger sorted file
usage:
   chainMergeSort file(s)
Output goes to standard output
options:
   -saveId - keep the existing chain ids.
   -inputList=somefile - somefile contains list of input chain files.
   -tempDir=somedir/ - somedir has space for temporary sorting data, default ./

liftover

liftOver -minMatch=0.9 -genePred Bos_taurus.ARS-UCD1.2.110.chr.genePred  ARS-UCD1.2ToOar_rambouillet_v1.0.chain Oar_rambouillet_v1.0.chain.gene.map Oar_rambouillet_v1.0.chain.gene.unMapped #直接复制过来的脚本

./liftOver -minMatch=0.9 inputfile ARS-UCD1.2ToOar_rambouillet_v1.0.chain outputfile unmappedfile
关于liftover的使用教程有很多，随便搜搜都是，最大的难点是chain文件