litover只能提供一些熱門基因組的轉換，而西農(nóng)開發(fā)的RGD只有一小部分基因組爆哑，不過非常棒的是他們提供了maf文件什猖，這是手動創(chuàng)建chain文件中最耗時間和資源的一步。本文主要記錄如何以西農(nóng)提供的文件獲取chain文件并最終做到坐標的轉換素邪。
當然如果RGD沒有提供所需的chain文件外莲，也可以參照他們的腳本自主創(chuàng)建，不過這一步需要耗費大量的時間和服務器資源（sheep1比對cattle1花費了我 24個小時和近200g的內存）

文件下載和準備

1 maf文件

http://animal.omics.pro/code/index.php/RGD
非常棒的數(shù)據(jù)庫兔朦，有很多有用的工具和文件

maf文件下載：RGD顯示 (omics.pro)
下載完畢后建議用md5校驗一下：

md5sum -c md5.txt 

該論文的附表中提供了上述maf文件的參考基因組信息：RGD v2.0：反芻動物功能和進化基因組學數(shù)據(jù)庫的重大更新 |核酸研究 |牛津學術 (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.常染色體和性染色體是數(shù)字形式摆舟，需要進行染色體號的轉換

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

后續(xù)又發(fā)現(xiàn)一個問題亥曹，maf文件中的染色體號并不是標準的數(shù)字染色體號，而是 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這個由西農(nóng)自主組裝的基因組）
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

# 檢查參數(shù)是否正確
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開頭的行，開始一個新區(qū)塊
        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
        # 判斷區(qū)塊中是否有以s Sheep開頭的行
        if [[ "$block" == *"s Sheep"* ]]; then
            # 輸出區(qū)塊的前兩行和以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#這里整個循環(huán)挺不錯的

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' 讀取數(shù)據(jù): 沒有那個文件或目錄）
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參數(shù)是必須的
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 #不懂參數(shù)，直接復制來的

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文件