上一次小编介绍了如何下载物种的contig数据，没看过的小伙伴可以点击以下链接：
https://www.jianshu.com/p/08104688516f
这次小编介绍下如何做Alignment以及SV 谱系分配（lineage assignment）？
第一步: 配置所需程序路径，配置前请确保你的机器上安装了以下软件：
hdf5，bzip，lzmalib，openssl，snakemake1，python，blasr，smartie-sv-master
配置过程：
export HDF5INCLUDEDIR=/path /hdf5/include
export HDF5LIBDIR=/path/hdf5/lib
export LD_LIBRARY_PATH=/path/hdf5/lib:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH=/path/zlib-1.2.7:/path/bzip2-1.0.6:/path/lzmalib-0.0.1:/path/curl-7.28.1/lib:/path/openssl-1.0.1:$LD_LIBRARY_PATH
export PATH=/path/snakemake/bin:$PATH
export PYTHONPATH=/path/python3.6/lib/python3.6/site-packages:/path/snakemake
第二步：构建基因组索引，在此我们用人类(GRCh38) & 黑猩猩（pantro6）为例
/smartie-sv-master/bin/sawriter panTro6.contig.fasta
/smartie-sv-master/bin/sawriterGCF_000001405.38_GRCh38.p12_genomic.fa
第三步：运行snakefile
# target：hg38 query：pantro6
snakemake -s Snakefile -w 50 -p -k -j 20
Snakefile：
shell.prefix("source config.sh; set-eo pipefail ; ")
configfile: "config.json"
def _get_target_files(wildcards):
return config["targets"][wildcards.target]
def _get_query_files(wildcards):
return config["queries"][wildcards.query]
rule dummy:
input: expand("variants/{target}-{query}.svs.bed",target=config["targets"], query=config["queries"])
rule callSVs:
message: "Calling SVs"
input :SAM="mappings/{target}-{query}-aligned.sam",TARGET=_get_target_files, PG=config["install"] +"/bin/printgaps"
output : "variants/{target}-{query}.svs.bed"
shell : """
cat {input.SAM} | {input.PG} {input.TARGET}variants/{wildcards.target}-{wildcards.query}
"""
rule runBlasr:
message: "Aligning query to target"
input: BL=config["install"] + "/bin/blasr",TARGET=_get_target_files, QUERY=_get_query_files
output: "mappings/{target}-{query}-aligned.sam","unmappings/{target}-{query}-unaligned.fasta"
shell: """
{input.BL} -clipping hard-alignContigs -sam -minMapQV 30 -nproc 6 -minPctIdentity 50 -unaligned{output[1]} {input.QUERY} {input.TARGET} -out {output[0]}
"""
config.json
{
"install" :"/path /smartie-sv-master",
"targets" : {
"hg38" : "/hg38/GCA_000001405.27_GRCh38.p12_genomic.retain.fasta"
},
"queries" : {
"pantro6" : "/Path/pantro6/panTro6.contig.fasta "
},
}
config.sh
# this is only if you don't have hdflibrary globally installed
exportLD_LIBRARY_PATH=$LD_LIBRARY_PATH:/net/eichler/vol5/home/mchaisso/software/hdf/lib
以上是如何进行Alignment以及Calling SV，下面介绍下如何做SV assignment？
首先介绍下Assignment原理
表格中1表示在物种中出现，0表示不出现，1-5表示物种1-5
如果SV1在物种1中出现，且只与物种2中某一个SV（任何一个都可以）重叠，则认为该SV是1,2shared SV。
如果SV1在物种1中出现，且与物种2和3中某一个SV（任何一个都可以）重叠，则认为该SV是1,2,3 shared SV。
表格中的其他lineage分配类似。
注意：
1）后期在分配时只考虑SV有没有重叠，不需要考虑他们的坐标，也不需要合并他们的坐标。
2）后期判断某一个SV在物种中有没有，就看他有没有和物种里的SV重叠，只要有一个重叠，就认为它在该物种中有
1, AssignDeletion缺失
第一步：得到每两两物种之间overlap的SV列表，
只要两个SV重叠度大于50%则认为是同一个SV（这里采用的标准是参照reference【1】里的标准，也可以设置其他的标准）
for i in `ls *del.bed`
do
for j in `ls *del.bed`
do
if[ $i != $j ]
then
bedtools intersect -a $i -b $j -f 0.5 -r -wo >${i%.*}.${j%.*}.common
fi
done
done
第二步：将SV分配到各个分支上去
perl convert5SpeciesSVToLineageSV-v2.pl--allSVList allSVList --commonSVFileMat commonMat-5 --outPrefix out > log
2. Assign insertion插入
首先你需要将SV_end =SV_start + svLen. 生成一个新的SV bed文件，然后参照缺失deletion的做法进行分配。
注：convert5SpeciesSVToLineageSV-v2.pl是小编自写脚本。如有需要，请进QQ群下载，群号：756263353， 文件分享在群文件里面
参考资料：
【1】2018-Science- High-resolution comparative analysis of great apegenomes
