趋同检测1
2023-11-12 本文已影响0人
多啦A梦的时光机_648d
总体来说分为四大步
##1. 利用rbh的方法是别ortholog
##2. 利用ortholog(每一个ortholog里面只保留物种名)和convgent_1.2.py生成prank文件和01_NodeSeq，02_Trees，03_rate，04_siteFreq文件，其中prank用于后续的pcoc输入，而01_NodeSeq，02_Trees，03_rate，04_siteFreq用于后续的convCal的输入
##3. 利用probCal.py脚本检测趋同信号
## 4. 运行pcoc.1.0.py脚本检测趋同
切记树文件必须是种属名这种，不能像正选择那样简写,这是为了后续pcoc可以直接进行
##脚本在https://github.com/BIGtigr或者https://github.com/taotaoyuan/
1. 利用rbh的方法鉴定ortholog
cd /home/lx_sky6/yt/0729_Carex/20-Conv
export PATH=/home/lx_sky6/software/miniconda3/envs/yt/bin:$PATH
##1-renamed_pep为一个文件夹，里面包含了需要比对的物种蛋白序列,名字需要改为*.aa.fas
python batch_reciprocal_blast.py REF.fas 1-renamed_pep
export PATH=/home/lx_sky6/software/miniconda3/envs/yt/bin:$PATH
python2 /home/lx_sky6/yt/0310_aquatic_botany/ZDH/parsingLIST.py ./ /home/lx_sky6/yt/0729_Carex/20-Conv/2-cds
##ortholog目录一共生成2212个ortholog.fasta
2. 利用convgent_1.2.py生成prank,convert_paml,run_paml和01_NodeSeq，02_Trees，03_rate，04_siteFreq文件，其中prank用于后续的pcoc输入，而01_NodeSeq，02_Trees，03_rate，04_siteFreq用于后续的convCal的输入
cd /home/lx_sky6/yt/0729_Carex/20-Conv/5-Converg
cp -r ortholog ortholog_new
sed -i 's/-.*//g' ./ortholog_new/*.fasta
```{ortholog22136.fasta}长这样
>Sorghum_bicolor
MQTAAAVSASPPDLLRRRLHRPAAAGLPSLSFSSSPSRNLGCLVSNTGSRFSRVVVASSPASFDGLDAKDSSISSSSNSSVLVSLVQEIEPLDVSLIQKEVPAETVDAMKRTISGMLGLLPSDQFKVLVEAPWDPVKRLLVSSMMTGYTLRNAEYRLSFERNLDFSEEDSNNKISRESAQTEASRSDELAPGSEKKDNNEENNKLNQDLGFEGLGELSSQARDYIIQMQARLNSIKKELHNLKRKNSALQMQQFVGEEKNDLLDYLRSLKPEKVAELSEPTCEGVKEAIHSVAHGLLATLSPKMHSIKNNNINNNDNNNAALNFGKNENNNDDDCTELVENTSINFQPHLSVPRDYLARLLFWCMLLGHYLRGLEYRLELMQLYNISKGLGSTPNNNNNGDFVA
>Zea_mays
MPTASSLRPPLPSSAARRLRFPPPGQLLPAGLAIGSRRRRLAGVGVAAASASPFDELYARGRPVHGPSKKSILWNLIQDIEPLDLTVIQKDVPPETVDAMKRTISGMLGLLPSDQFRVVVEALWNPFFKLLVSSIMTGYTLHNAEYRLSFERNLELSEEDAEFQRRDITEDNHNDINLGRPVTIFRLSEEGMSQDSGKSDEESCCESMGEELGNLTPQAEEHIIRLQSRLDAMEKELHDLKRKNSALQMQQFVGEEKNDLLDYLRSLTPEKVAELSESTCPGVQEAIQSVVHGLLATLSPKIHSKSPPPLENAAGGALNLGGENDDCAELVENASLPFQPLISVPRDYLARLLFWCMLLGHYIRGLEYRLELAQLLRISSDVGSFSGGNDHVI
运行convgent_1.2.py
11 1
(((((Sorghum_bicolor,Zea_mays),((Brachypodium_distachyon,Achnatherum_splendens),Oryza_sativa)),((Juncus_inflexus,Juncus_effusus),(Carex_breviculmis,Carex_littledalei))),Ananas_comosus),Arabidopsis_thaliana);
python convgent_1.2.py --path=/home/lx_sky6/yt/0720_GJ_genome/00.Converg/4-Converg/ortholog_new --tree=/home/lx_sky6/yt/0720_GJ_genome/00.Converg/4-Converg/tree.tre
##很奇怪，不知道为什么这里只会出来prank，convert_paml，run_paml结果，后面的结果不能一次出来，所以需要注释convgent_1.2.py中前面的os.system()，然后继续进行后续的，才会生成01_NodeSeq，02_Trees，03_rate，04_siteFreq
##生成文件
prank
convert_paml
run_paml
01_NodeSeq
02_Trees
03_rate
04_siteFreq
@@@@@@@@注意@@@@@@@@@@#########
注意1: 检测01_NodeSeq和02_Trees个文件中是不是有0内容的，要删除，不然会报错（虽然genelist.txt是已经去掉了为0的基因了，但是为了稳妥自己再找一次）
cd 01_NodeSeq
find . -type f -empty|sed 's/.///g'|sed 's/.nodes//g' >../empty.list ##241个
##-rw-rw-r-- 1 lx_sky6 lx_sky6 0 Oct 13 00:40 ortholog10796.phy.ntree
##-rw-rw-r-- 1 lx_sky6 lx_sky6 0 Oct 13 00:40 ortholog10602.phy.ntree
在最后要计算的list中去除即可
注意2: 由于序列差异太大，会导致保守的氨基酸位点太少，paml无法计算rate，会生成-nan，导致最后一步报错，这些-nan的基因就从testGeneList中删除。
cd 03_rate
for id in $(cat ../genelist.txt); do grep 'nan' $id.rat |awk -v var="$id" '{print var}' ;done
会得到几个结果为nan的基因名，记得把这几个和empty.list都存入remove.id中删除 ##这次没有
################@@@@@@@@@@@##################
最后
grep -vf remove.id genelist.txt >testGeneList
