从参考基因组(序列)和SNPs(位点)入手的Ka/Ks（dN/d
网上查了查似乎没有一个完整的从参考基因组reference.fna和注释文件.gff出发，到blastn，再到提取orthologs，再到ka/ks选择分析的流程。这里记录一下自己的实验过程，理论上可行（不做任何担保）。
1、下载fna参考基因组以及gff注释文件
不详细说了，NCBI搜到以后wget就行。
2、根据.gff文件提取联合的CDS区域
这一步可以自己写脚本，也可以用gffread这个软件
github主页：https://github.com/gpertea/gffread
使用很简单
gffread reference.gff3 -g reference.fna -x cds.fa
关于CDS、exon、mRNA之间的区别我之前写过。简单来说sum(exon)-UTR=CDS, sum(exon)+ployA+G帽子=mRNA。
所以如果看.gff文件会发现每一个mRNA都对应包含很多exon，也对应包含与exon同样多的CDS，我们要的是编码蛋白的CDS。因为在Ka/Ks分析中，需要保持”是否编码蛋白“这一标准相同，编码区和非编码区差别的Ka/Ks差别很大。最后可以看到gffread提取后的CDS是合并到一起的，在.gff文件中它们的ID都相同。
3、blastn
先对要提取orthologs的每一个物种的CDS进行建库，然后两两比对，双向比对bit score都为最高的被认为是这两个物种的一对ortholog。但是如果物种太多的话两两比对工作量会爆炸性增长，所以可以确定一个物种作为比对对象，其他几个物种和它进行两两比对。比如有ABCDE五个物种，确定A为标标准准对象，进行A-B,B-A,A-C,C-A,A-D,D-A,A-E,E-A，10次比对。筛选过后再以A为对象进行merge合并。
先对5个物种的CDS进行建库，假设文件名储存在name.list文件里：
for i in $(<name.list);do makeblastdb -in ${i}.CDS.fasta \
-dbtype nucl -parse_seqids ;done
以其中之一(比如说PT物种)为标准进行双向的比对：
for i in $(<name.list);do { blastn -db PT.CDS.fasta -query ${i}.CDS.fasta \
-outfmt 7 -out ${i}_to_PT_blast ; }& done
for i in $(<name.list);do { blastn -db ${i}.CDS.fasta -query PT.CDS.fasta \
-outfmt 7 -out PT_to_${i}_blast; }& done
并行十分钟差不多能跑完
4、初步处理blastn输出文件
这里就仁者见仁智者见智了，例如
cat XXX_to_PT_blast.out|grep -v '#'>XXX_to_PT_blast.out.edited
## 去除注释信息
for i in `cat XXX_to_PT_blast.out.edited|awk '{print $1}'|uniq|head`;do \
cat NL_to_PT_blast.out.edited|awk '$1=='"${i}"'{print $0}'|sort -k12,12nr| \
head -3>>${i}.best3;done
## 对每个query序列的blast结果筛选最佳的(bit score)三次比对
这样就筛掉了很多分数较低的比对，也去掉了部分旁系同源基因paralogs的干扰。后续还需要进一步筛选。这里没有选择唯一最高bits core基因，而是使用了最高三个的原因是，一些CDS的bit score是相同的，还不能分辨谁是最好的，要看另一个方向的blast。
5、写个python脚本把两个方向blast的结果整合一下，选取双向配对的orthologs对
from multiprocessing import Pool
namelist=["BB","CC","DD","EE"]
def deal(name):
name=name.strip()
with open("PT_and_%s_orthologs"%name,"a+") as writefile, \
open("PT_to_%s_blast.best3"%name,"r") as foward, \
open("%s_to_PT_blast.best3"%name,"r") as backward:
fowardlist=foward.readlines()
backwardlist=backward.readlines()
for line1 in fowardlist:
line1_spl=line1.strip().split("\t")
for line2 in backwardlist:
line2_spl=line2.strip().split("\t")
if line1_spl[0]==line2_spl[1] and line1_spl[1]==line2_spl[0]:
## 如果某次方向1的第一列(query)等于某次方向2的第二列(database)，并且
## 方向2同一次blast的第一列(query)存在于方向1同一次blast的第二列(database)，
## 认为其为最佳匹配
writefile.write(line2)
continue
if __name__ == '__main__':
pool = Pool()
pool.map(deal, namelist)
pool.close()
pool.join()
## 多线程
6、选取唯一的orthologs对
现在我们得到的orthologs对仍然是不唯一的，因为我们之前使用的是bit score top3的筛选策略，难免出现paralogs也在top3 范围内双向匹配成功的情况。这个时候，第一是根据bit score降序，第二是根据alignment长度降序排列，选出top1来。
cat PT_and_XXX_orthologs|awk '{print $1}'|uniq>PT.CDS.name.list
cat PT_and_XXX_orthologs|awk '{print $2}'|uniq>XXX.CDS.name.list
for i in $(<PT.CDS.name.list);do cat PT_and_XXX_orthologs| \
awk '$1=''"${i}"''{print $0}'|sort -k12,12 -k4,4|head -1>> \
PT_and_XXX_orthologs.1;done
## 第一个方向的uniq
for i in $(<XXX.CDS.name.list);do cat PT_and_XXX_orthologs| \
awk '$2=''"${i}"''{print $0}'|sort -k12,12 -k4,4|head -1>> \
PT_and_XXX_orthologs.uniq;done
## 第二个方向的uniq
7、5个物种merge
现在每个PT_and_XXX_orthologs.uniq文件中的CDS关系都是一对一的了，即一个来自PT的CDS只有一个来自XXX的CDS与之匹配，反之依然。
因为我们是将PT这个物种作为参考，所以merge的时候也用PT物种。
我这里用的是R脚本
data1<-read.table("BB_and_PT_orthologs.uniq")
data2<-read.table("CC_and_PT_orthologs.uniq")
data3<-read.table("DD_and_PT_orthologs.uniq")
data4<-read.table("EE_and_PT_orthologs.uniq")
data1<-data1[,1:2]
data2<-data2[,1:2]
data3<-data3[,1:2]
data4<-data4[,1:2]
title1<-c("AA","PT")
title2<-c("BB","PT")
title3<-c("CC","PT")
title4<-c("DD","PT")
colnames(data1)<-title1
colnames(data2)<-title2
colnames(data3)<-title3
colnames(data4)<-title4
data_merge1<-merge(data1,data2,by="PT",all=FALSE)
data_merge2<-merge(data_merge1,data3,by="PT",all=FALSE)
data_merge3<-merge(data_merge2,data4,by="PT",all=FALSE)
table(duplicated(data_merge4[,1]))
table(duplicated(data_merge4[,2]))
table(duplicated(data_merge4[,3]))
table(duplicated(data_merge4[,4]))
table(duplicated(data_merge4[,5]))
data_merge4<-data_merge3[-1913,]
## 检查一下有没有重复的，理论上应该没有，但可能也会出现一两个，可能是计算上出了小误差。把重复的删掉就行。
write.table(data_merge4,"all_5_orthologs",quote=F,row.names = FALSE)
8、从SNPs到orthologs
比如说现在有一个物种KK的SNPs，已经被call出来了，而且参考基因组是BB。现在想把它一起放进去做选择分析的话可以先把SNPs map回reference（这个reference一定要是5个物种之一，在本例中为BB），得到alternate reference fasta，然后根据位置的对应关系，提出orthologs。
第一步，先进行筛选，超过群体数量/2的个体占有的纯合SNPs(也就是1/1)才被考虑
import pandas as pd
from sklearn import datasets
my_data = pd.read_table('BB.orth.seq.CDS.loc')
df_orth=pd.DataFrame(my_data)
my_data2=pd.read_table("filter.mod.HN.vcf")
df_SNPs=pd.DataFrame(my_data2)
if_hetero_proper=[]
for i in range(0,len(df_SNPs)):
num=0
for k in range(9,19):
if df_SNPs.iloc[i,k].startswith("1/1"):
num=num+1
elif df_SNPs.iloc[i,k].startswith("0/0"):
num=num+0
else:
num=num-99
if_hetero_proper.append(num)
df_SNPs['hetero_flag']=if_hetero_proper
df_SNPs.to_csv('filter.mod.HN.flagged.vcf',sep='\t', index=False)
上面这步是加上flag。下一步根据flag筛选就行。
cat filter.mod.HN.flagged.vcf|awk '$20>=5||$20=="hetero_flag" \
{$20=NULL;print $0}'>filter.mod.HN.flagged.filtered.vcf
因为要对project保密，有一些文件名会稍显混乱，大家见谅，交流思路。
下一步就把筛选出的SNPs mapped回去。可以用GATK的FastaAlternateReferenceMaker
但是在此之前先把筛选过的VCF处理一下，把空格分隔符改为tab分隔符，再加个VCF版本号，再index一下
cat filter.mod.YN.flagged.filtered.vcf|sed 's/ /\t/g'|sed '1i\ ##fileformat=VCFv4.2'> \
filter.mod.YN.flagged.filtered.tabbed.vcf
gatk IndexFeatureFile -I filter.mod.YN.flagged.filtered.tabbed.vcf
remap
gatk FastaAlternateReferenceMaker -R BB.fasta \
-O KK.SNPs.remapped.fasta -V filter.mod.YN.flagged.filtered.tabbed.vcf
## gatk FastaAlternateReferenceMaker \
## -R reference.fasta \
## -O output.fasta \
## -L input.intervals \
## 如果需要在固定区间上map回去可以选择这个参数，非必要
## -V input.vcf \
FastaAlternateReferenceMaker得到的序列索引有点问题，比如说本来应该长这样的>NC_01010101.1 ，它现在长这样>2 NC_01010101.1:1-123690457。需要改回来
cat YN.SNPs.remapped.fasta｜grep '>' >reads.name
for i in $(<reads.name);do k=`grep ''"${i}"'' KK.SNPs.remapped.fasta`; \
echo ${k};sed -i 's/'"${k}"'/>'"${i}"'/g' KK.SNPs.remapped.fasta;done
9、根据 orthologs（CDS） ID 信息提取序列
## 对remapped的物种KK重复第二步，提取出CDS序列
gffread BB.gff -g kk.fasta -x KK.CDS.fasta
awk 'NR==FNR{ids[$0];next} /^>/{f=($1 in ids)} f' BB.orthologs.id \
KK.CDS.fasta > KK.orthologs.uniq.fasta
## 这个是KK物种的remapped SNPs的orthologs，因为其call SNP的时候就是以BB为参考基因组的，\
## 所以直接把BB的orthologs序列ID用来提KK的序列就可以了。
## 同理将所有其他物种的序列提出来，依据其ID。ID我们之前通过第7步已经得到了。
for i in "AA" "BB" "CC" "DD" "EE" "PT"
do
{
awk 'NR==FNR{ids[$0];next} /^>/{f=($1 in ids)} f' ${i}.orthologs.id \
${i}.CDS.fasta > ${i}.orthologs.uniq.fasta
}&
done
10、orthologs alignment
假设现在是6个物种（AA,BB,CC,DD,PT,KK），6个文件，每个文件一万多条fasta序列。因为要做alignment，我们将其转换为一万多个文件，每个文件6条序列（按顺序）。
for i in $(<name.list)
## name.list 是包含六个物种名称的文件，也就是（AA,BB,CC,DD,EE,KK）
do
{
for((count=1;count<=13435;count++))
do
{
name=`cat ${i}.name.list|sed -n ''"${count}"'p'`
## ${i}.name.list 是每个物种的CDS ID文件
cat ${i}.orthologs.uniq.fasta|seqkit grep -p ''"${name}"'' \
>>./sep_orthologs/combine_orthologs.${count}
## 需要安装seqkit
}
done
}&
sleep 30
done
一些其他可能用到的操作
for i in {1..13435};do awk '/^>/{print ">" ++i; next}{print}' combine_orthologs.${I}> \
../sep_orthologs_namechanged/combine_orthologs.namechanged.${i};echo ${i};done
## 把抬头改为>1 >2等等
seqkit concat --quiet combine_orthologs.namechanged.*> \
../concatenated.orthologs.fasta
## 串联所有抬头相同的序列(concatenated)
alignment这里可能是个坑。因为Ka/Ks分析以及codeml的基本原理都是基于codon（三联体密码子）的，所以要求序列为3的倍数。常用的clustalo似乎没有这个功能。在线版的MAFFT好像可以，但是本地版的我没找到这个功能。最后使用了2018年发布的MACSE。它是基于codon进行alignment的，现生成AA序列，然后alignment，然后对应回核酸的alignment序列，而且会考虑终止密码子等等。比较适合Codeml。（注意:如果你有很长的序列，一定要记得设置java参数 -Xmx，比如100G内存 -Xmx100G，否则会终止。）
macse：https://academic.oup.com/mbe/article/35/10/2582/5079334
使用macse对每个orthologs进行codon-alignment。
for k in $(<ortho_sep09)
do
{
java -jar \
~/software/macse.jar -prog alignSequences \
-seq ./sep_orthologs/${k} \
-out_NT ./macse_codon_align/${k}
}
用代码检验一下是不是所有orthologs位点数都为3的倍数：
for i in {1..13435};do k=`cat codon_align_combine_orthologs.namechanged.${i}| \
sed -n '2p' |sed -e 's/-//g'|wc -L`;echo `expr "scale=2;${k}/3"|bc`;done
得到结果：
224.00
405.00
1151.00
625.00
293.00
582.00
...
都能整除3，到位了。说明这个软件可行。
记录一下存在移码突变的orthologs（是用!标注的位点），以后可能要排除：
for i in {1..13435};do k=`cat codon_align_combine_orthologs.namechanged.${i} \
|grep '!'`; p=$((p+1));if [[ -n $k ]];then echo $p>> \
../frameshift.list; fi; done
11、 Ka/Ks （dN/dS）选择分析
使用的还是最nb(也最笨重)的杨子恒老师的PAML包，里面的codeml
PAML包非常吃数据格式。它像是一个笨重而缓慢的巨人，一个穿戴满金银首饰、令人眼花缭乱的贵人，也像是一个有无数支线的副本。所以如果我想得到每一个orthologs的Ka/Ks应该怎么处理数据呢？如下
orthologs输入格式示范
在这个文件中，有6个物种，3个(组)orthologs同源基因，这三个orthologs-alignment的长度分别为20，40，60（当然了这个数字是假的，真实输入的一定要是3的倍数！），六个物种的名称分别为AA,BB,CC,DD,PT,KK（真实的名称应该不能取ATCG这种，否则可能读取错误）。每一个alignment其实都是phylip格式的，把它们cat到一起输入就行了。
也可以是这样的
也可以是这样的，这里1-7是名称，空两格以上接序列，不用换行。这里显示的只是一个orthologs。
可能需要用到fasta alignment转换为phylip alignment的工具，这里有一个perl脚本：https://github.com/liaochenlanruo/catfasta2phyml
注意避坑：“一些错误”
记得去除移码突变的loci
codeml配置文件
seqfile = all.final.phy * sequence data filename
treefile = codeml.label.HN.tree * tree structure file name
outfile = mlc_HN * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 0: concise; 1: detailed, 2: too much
runmode = 0 * 0: user tree; 1: semi-automatic; 2: automatic
* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise
seqtype = 1 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
ndata = 13435
clock = 0 * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
aaRatefile = dat/jones.dat * only used for aa seqs with model=empirical(_F)
* dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own
model = 2
* models for codons:
* 0:one, 1:b, 2:2 or more dN/dS ratios for branches
* models for AAs or codon-translated AAs:
* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
* 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)
NSsites = 2 * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
* 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
* 13:3normal>0
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
Mgene = 0
* codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
* AA: 0:rates, 1:separate
fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = .4 * initial or fixed omega, for codons or codon-based AAs
fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
Malpha = 0 * different alphas for genes
ncatG = 8 * # of categories in dG of NSsites models
getSE = 0 * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 1 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)
Small_Diff = .5e-6
cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?
* fix_blength = 1 * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
method = 0 * Optimization method 0: simultaneous; 1: one branch a time
* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK.
建议做两个模型，一个选择模型，一个零假设模型（neutral evolution）。这里使用branch-site模型A，也就是model=2，NSsite=2, fix_omega=0
具体选项复杂，一些是可以默认的，还是必须得看说明书。说明书里写了这个模型比较适合选择分析。
然后处理得到的模型似然值，计算零假设和备则假设之间的2delta(lnL)，带入自由度为1的卡方分布或50:50原点和卡方混合分布，得到p值，再根据p值进行Benjamini-Hochberg算法矫正，得到q值。设定q阈值0.1, 0.01, 0.05得到受选择基因。
branch-site model A
branch-site model A
树文件，建议用raxml跑：
(((((NL,YN),HN #1),(PT,PA)),MM),CJ);
运行：
nohup codeml codeml.ctl &
2021.1.6 更新
发现其实提取orthologs不用这么麻烦，直接用last做两两比对就行了。
last太强了。
web: http://last.cbrc.jp
tutorial1: http://last.cbrc.jp/doc/last-tutorial.html
tutorial2: https://github.com/mcfrith/last-genome-alignments
