（13）共线性比对相关（感谢zyr师姐大大！！！）

#输入文件：

1）待比对的全基因组文件(屏蔽了重复序列的)

2）参考物种的gff格式文件（gene transfer format）,用来对基因组进行注释

#一定要注意！基因组文件里面的scaffold和gff3的命名要一致，这样才能切基因，做下游！！！
####再次注意！！！做基因预测的时候字符不能大于16，所以还是拿scaffold_XX，去预测，然后注意基因组文件和gff3的对应！！！

后期有空再出师姐的代码解读工作

一、准备工作

mkdir 0828bidui #在~/Qxy/knowngenome/gongxianxing目录下
首先需要保证24个基因组＋1个外群都已经做了重复序列分析，有Repeat_result，然后对其中的软屏蔽基因组文件进行改名，基因组序列前面加上物种名，然后后缀都改成*_genomic.fna统一格式

返回到脚本目录下修改脚本
vim changename.sh #在~/Qxy/qxyjiaoben文件夹下更改shell脚本，改对了用于更改文件中序列名

返回到repeat_result文件夹下去修改基因组序列名字
sh ~/Qxy/qxyjiaoben/changename.sh
mv *_genomic.fna ~/Qxy/knowngenome/gongxianxing/0828bidui/

二、开始运行比对切得基因拿到四倍简并位点

2.1跑lastZ

补充比对小知识：一般来说，重测序数据比对多采用BWA和Bowtie2;
不同物种间基因组（不同物种的参考基因组）的比对常用共线性比对，如LastZ，此外还有Last

python3 querysplit_test.py ref_genomic.fna query_genomic.fna absolute_path_of_genom_file_directory(/ifs1/User/wuqi/...)
输入参数：1）python脚本 2）参考基因组 3）待比对基因组 4）工作目录的绝对路径
！！！注意啊大姐！！！这里输入的参数是要做共线性比对的参考基因组，而不是所有基因组，换言之就是想让其他的都跟XX比，就选XX做参考基因组 so南极石耳是参考基因组，其他的是待比对的，包括外群！

!!!parafly批量运行!!!
vim parafly.txt
文件内每一行写入这样一条命令 
python3 /ifs1/User/wuqi/Qxy/qxyjiaoben/fromzyr/querysplit_test.py Umbilicaria_antarctica_genomic.fna Hypocenomyce_scalaris_genomic.fna /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0828bidui
就是将南极石耳与其他24个（含外群）进行比对
nohup ParaFly -c parafly.txt -CPU 40 & 
挂后台批量运行，25个大概两天左右
结果文件是参考基因组有.2bit .size .spil 其他基因组有 .2bit .size .split 多了个.axt
这个axt目录下就存放着多个以该基因组的scaffold序列命名的文件夹，再进去，里面就有maf结果文件

2.2跑multiz

补充知识：

在某些情况下，使用 Lastz 进行初步比对，然后再使用 Multiz 进行进一步的多序列比对可以提供更全面和准确的结果。

Lastz 是一个广泛应用于序列比对的工具，它具有高速和较高的灵敏度。通过使用 Lastz，可以在两个序列之间执行较快的局部比对，找到相似区域和保守区域。这种初步比对可以帮助确定序列之间的相似性和重要的结构/功能区域。然而，Lastz 主要针对两个序列的比对，并且在涉及多个序列同时进行比对时可能不够高效。这时候就可以引入 Multiz。

Multiz 是专门设计用于多序列比对的工具，它能够同时比对多个序列，包括多个物种或同一物种内的多个基因组序列。Multiz 使用其中一条序列作为参考，在该参考序列的背景下将其他序列逐个与之比对。通过这种方式，Multiz 会生成一个多序列比对结果，展示了各序列之间的相对关系、共同保守区域等信息。

因此，首先使用 Lastz 对两个序列进行局部比对，可以快速发现重要的相似区域。然后，将这些序列与其他相关序列一起放入 Multiz 进行多序列比对，可以获得更全面的结果，包括多个物种或同一物种内的多个序列之间的共同保守区域，进化关系等。这样的策略可以综合利用两个工具的优点，提供更全面和准确的多序列比对分析。

python3 newmultiz.py ref_genomic.fna workdictory(共线性比对目录) 

nohup python3 /ifs1/User/wuqi/Qxy/qxyjiaoben/fromzyr/newmultiz.py Umbilicaria_antarctica_genomic.fna /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0828bidui &
用newmultiz.py脚本，输入参数是参考基因组文件名和基因组文件所在目录
11个文件要跑三个小时左右，25个文件大概跑了一天一夜

#这个运行完会生成chr_multiz/文件夹，在该目录下有以参考基因组的序列为名的文件夹，每条序列一个，里面会含有maf文件

2.3跑切基因

python3 /ifs1/User/wuqi/zyr/convergent/maf2gene_newum.py ref_genome.fna path_of_genom_file_directory(/ifs1/User/wuqi/...) absolute_path_of_gff

#运行maf2gene_newum.py这个脚本，输入参数是 1）参考基因组名称 2）本次进行共线性比对的基因组所在目录绝对路径 3）参考基因组的gff3文件绝对路径

nohup python3 /ifs1/User/wuqi/Qxy/qxyjiaoben/fromzyr/maf2gene_newum.py Umbilicaria_antarctica_genomic.fna /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0828bidui /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/Umbilicaria_antarctica_genomic.gff3 &

一定要保证XX.fna和XX.gff3文件中序列名字是一致的，这样才能切基因下来！！！
本次南极石耳预测到8391个基因，运行完毕后得到newgenes目录，该目录下以U.antarctica的每条scaffolds序列为名生成对应的目录，然后在该目录下有的被切出来好几条基因，有的则没有被切出来
25个大概跑了半天左右

注！！！这次石耳切基因出现问题，是由于每行末尾不是；代码就报错了，以后注意！！！

2.4然后过滤一下

python3 ~/zyr/convergent/filter.py refgenome workdirectory

#运行filter.py这个脚本，输入参数是参考基因组文件和此次共线性比对工作目录

python3 /ifs1/User/wuqi/zyr/convergent/filter.py 工作目录 外群基因组序列文件 比对不上的物种数
#该代码运行完毕后会生成gapfiltered文件夹，其中是多个基因文件

python3 filter.py Umbilicaria_muehlenbergii_genomic.fna /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0926bidui Umbilicaria_antarctica_genomic.fna Hypocenomyce_scalaris_genomic.fna 0 

ref = sys.argv[1]
workdir = sys.argv[2]
foreground = sys.argv[3]#species name under branchsite model or none under site model
og = sys.argv[4]
allowednum = sys.argv[5]

##注意，这个脚本有更新！！！现在是五个参数1）参考基因组名称 2）工作目录 3）前景枝 4）外群 5）允许没有比对上的物种数（含外群）

如：我有24+1个物种，想得到含外群有3个物种没有比对上的，就可以设置最后一个参数为3 ，这样拿到的是剩下22个物种比对结果

已经对师姐脚本都创建了软连接~~~既不怕丢，也可以同步更新了

nohup python3 filter.py Umbilicaria_muehlenbergii_genomic.fna /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0926bidui Umbilicaria_antarctica_genomic.fna Hypocenomyce_scalaris_genomic.fna 0 &
#设置0 就表示不允许有哪个物种没有比对上，所以是每个物种都有的单拷贝基因
#设置其他数字，那就是有没比对上的

2.5获取四倍简并位点

四倍简并位点（Fourfold Degenerate Synonymous Site, 4DTv）,在进化学上被作为评估基因组是否发生全基因组复制事件的参数。

python3 4dtvseq.py *_genomic.fna的目录 含gapfiltered的目录

#运用4dtvseq.py脚本获取四倍简并位点序列，输入参数两个1）_genomic.fna文件所在的目录 2）gapfiltered文件夹所在目录

python3 /ifs1/User/wuqi/Qxy/qxyjiaoben/fromzyr/4dtvseq.py /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0828bidui /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0828bidui
这里改了一下师姐的代码，在sps.remove(ref.split('.')[0])前面加了一行if ref.split('.')[0] in sps:  这样就不会报错了
很快，就得到一个4dtv.fa文件，里面每个物种一条序列，其中第一条序列是我们的参考基因组序列，其名字不对，手动改一下，以物种名命名

三、建树建树建树！！！

3.1mafft比对（优先选择对话框输入）

mafft --auto input.fasta > output.fasta

#之前orthofinder建树时比对用的是muscle，速度快，但是只能比对短的序列（长度<500,序列数<500）

#这次muscle比对出来是空，所以选择mafft来比对

nohup mafft --auto --thread 20 4dtv.fasta > aln.4dtv.fasta &

#需要将XX.fa改为XX.fasta格式！！！
##补充##
服务器对话框直接输入mafft，可以进入到互动操作界面，依次选择以下5项，生成代码，跑起来
1）输入文件的路径，需要标准的fasta格式 如：/ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/4dtv.fasta
2）输出文件的路径 /ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/aln4dtv.fasta
3）选择输出文件的格式，输入对应序号按回车即可，支持Clustal、Fasta和Phylip三种格式。选3或视情况而定
4）选择比对策略，输入对应序号按回车即可，通常选择1
5）附加参数，如果没有，直接按回车即可
做完以上5个步骤，MAFFT工具通过你的各个选项，自动生成了运行命令，此时输入"Y"即可开始进行比对（其实该操作只是为了引导你生成代码）
#注，本次25个序列比对运行时间为12个小时左右，比较久，而且mafft跑命令行的时候结果文件前面出现了很多行别的东西，遂用了mafft对话框跑的内容

3.2Gblocks 去除空白序列

这个前面做一步将比对后的序列转换为单行
conda activate Gblocks #激活
Gblocks 4dtvbidui_no_blank.fasta -t gblocks.4dtv.fasta -b4=5 -b5=h -t=p -e=.2
这是Gblocks单命令行运行的代码，设置输入和输出文件就行。
输入文件是4dtvbidui_no_blank.fasta，输出文件是*.2和*.htm，跟gblocks.4dtv.fasta好像木得关系
这一步会生成*.2和*.htm的文件，而且会分段，cahtgpt说分段可能是为了保证对齐

3.3seqkit 去分段

seqkit seq 4dtvbidui_no_blank.fasta.2 -w 0 > gblocks.4dtv.fasta 
然后这一步会把.2文件转为单行，但是之间还是有空格
sed -i 's/ //g' gblocks.4dtv.fasta 
这一步是把空格去除掉

3.4raxml建树

#最终保留了4dtv.fasta、4dtvbidui_no_blank.fasta、gblocks.4dtv.fasta即原始4dtv数据，比对后的，比对后去空白序列的
#因为感觉空白序列很少，所以两者都建树，看看有区别不，一个服务器跑，一个在CIPRES上跑

nohup raxmlHPC-PTHREADS-SSE3 -f a -x 123456 -p 123456 -s 4dtvbidui_no_blank.fasta -m PROTGAMMALGX -N 1000 -n output -k -T 20 &

结果，服务器上建树超级超级超级慢，师姐是设-N 200都要跑一天，最终我24个物种跑了17天，结果文件如下：

1）程序执行信息已写入：RAxML_info.output

2）所有1000个引导树已写入：RAxML_bootstrap.output

3）得分最高的ML树已写入：RAxML_bestTree.output

4）带有支持值的最佳ML树已写入：RAxML_bipartitions.output 用的是这个树！！！

5）将支持值作为分支标签添加到最佳ML树中并写入：RAxML_bipartitionsBranchLabels.output

建树得到的结果和在线建树一致，就是U.caroliniana的枝系的支持值有点低，说是因为U.caroliniana是Umbilicaria和Lasallia的边缘种。