2023-03-21 | 渗透碎片检测 - Loter
渗入片段检测（单倍型）
基因组上的哪些片段是通过基因交流获得的以及它们的生物学功能
●适应性基因渗透包含两个过程：基因渗透和适应。
●基因渗透是指遗传成分从一个群体通过有性生殖的方式进入到另一个群体。产生杂交的前提是两个群体之间没有生殖隔离或只有不完全生殖隔离。只有这样才能通过与可育Ｆ１代回交达到遗传成份垂直传递的目的。
●通过基因渗透得到的区域，如果在该群体内没有受到选择，那么会因为随机遗传漂变的影响，逐渐在基因组中消失。只有当基因渗透得到的区域受到自然或人工选择，它才能在基因组中保留一定的频率。
软件对比：
与HAPMIX, LAM-LD, RFmix需要设置多个参数，HAPMIX需要提供遗传图谱、充足率、突变率和平均代时等较难获得的生物学参数。Loter除了单倍型数据以外不需要其他任何生物学参数即可完成祖先血统推断。与HAPMIX, LAM-LD, RFmix软件相比，在考虑人工模拟和人工混群的情况下，随着基因渗入次数的增加，Loter软件收到的影响较小。
Loter软件：
提供足够的祖先物种的基因组数据来作为来源群体（source populations），同时也需要发生渗入样本的基因组数据。对于每一个渗入个体的每一个变异位点，loter均会进行祖先来源判断，从所提供的来源群体中选择一个可能性最大的群体作为其来源。
需要提供一个几乎没有发生渗入或者渗入比例极低的群体作为对照。如将MSU即单倍型1视作对照，则结果文件中的1为未发生渗入区域
可以使用Dsuite的Fbranch模块的可视化结果得到
1、生成npy格式的文件
/home/sll/21-24.script/vcf2npy.py
python /home/sll/21-24.script/vcf2npy.py --vcffile x.vcf --samplelist sample.txt --outprefix x.npy
2、Loter推断单倍型渗入
loter_cli -r ref_phase.npy -a admix_phase.npy -f npy -o interval.txt
-r REF [REF ...], --ref REF [REF ...] files storing input reference haplotypes.可以是多个，空格分开输入
-a ADM, --adm ADM file storing input admixed haplotypes.
-o OUTPUT, --output OUTPUT
file to store the infered ancestries of admixed haplotypes, saved as a numpy array by default, or as a text CSV file if the file extension is '.txt'.
-f 输入文件的格式；有npy(默认), txt（csv）, vcf可选，所有的输入文件格式一致
-n CPU数量
-pc Run the phase correction module after the inference algorithm (only available for data from diploid organisms)
软件运行过程的日志：
loading file ANG.npy
Dimension of reference haplotype matrix 0 = (10, 52296)
loading file MSU.npy
Dimension of reference haplotype matrix 1 = (16, 52296)
loading file CHA.npy
Dimension of reference haplotype matrix 2 = (26, 52296)
说明，软件自动分渗入来源群体的单倍型，与最后的csv表格中的相对应
3、结果转换
/home/sll/21-24.script/npy2csv.py
python /home/sll/21-24.script/npy2csv.py --npyfile intrval.npy --outfile intrval.csv
两处python脚本均引自https://github.com/silvewheat/biotrans_cli，这是一位西农的师兄的github账号
vcf2npy.py
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 10 21:29:52 2020
@author: YudongCai
@Email: yudongcai216@gmail.com
"""
import click
import allel
import numpy as np
def vcf2npy(vcffile, samples):
callset = allel.read_vcf(vcffile, samples=samples)
haplotypes_1 = callset['calldata/GT'][:,:,0]
haplotypes_2 = callset['calldata/GT'][:,:,1]
m, n = haplotypes_1.shape
mat_haplo = np.empty((2*n, m))
mat_haplo[::2] = haplotypes_1.T
mat_haplo[1::2] = haplotypes_2.T
return mat_haplo.astype(np.uint8)
@click.command()
@click.option('--vcffile')
@click.option('--samplelist')
@click.option('--outprefix')
def main(vcffile, samplelist, outprefix):
"""
convert vcf to Loter input
"""
samples = [x.strip() for x in open(samplelist)]
mat_haplo = vcf2npy(vcffile, samples)
np.save(f'{outprefix}', mat_haplo)
if __name__ == '__main__':
main()
npy2csv.py
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 25 15:01:41 2020
@author: YudongCai
@Email: yudongcai216@gmail.com
"""
import click
import numpy as np
@click.command()
@click.option('--npyfile')
@click.option('--outfile')
def main(npyfile, outfile):
array = np.load(npyfile)
np.savetxt(outfile, array, delimiter='\t', fmt='%.0f')
if __name__ == '__main__':
main()