上一次小編介紹了如何下載物種的contig數(shù)據(jù)，沒看過的小伙伴可以點(diǎn)擊以下鏈接：
http://www.reibang.com/p/08104688516f
這次小編介紹下如何做Alignment以及SV 譜系分配（lineage assignment）震鹉？
第一步: 配置所需程序路徑俱笛，配置前請確保你的機(jī)器上安裝了以下軟件：
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
第二步：構(gòu)建基因組索引，在此我們用人類(GRCh38) & 黑猩猩（pantro6）為例
/smartie-sv-master/bin/sawriter panTro6.contig.fasta
/smartie-sv-master/bin/sawriterGCF_000001405.38_GRCh38.p12_genomic.fa
第三步：運(yùn)行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
以上是如何進(jìn)行Alignment以及Calling SV宽涌，下面介紹下如何做SV assignment平夜？
首先介紹下Assignment原理
表格中1表示在物種中出現(xiàn)，0表示不出現(xiàn)卸亮，1-5表示物種1-5
如果SV1在物種1中出現(xiàn)忽妒，且只與物種2中某一個SV（任何一個都可以）重疊，則認(rèn)為該SV是1,2shared SV兼贸。
如果SV1在物種1中出現(xiàn)段直，且與物種2和3中某一個SV（任何一個都可以）重疊，則認(rèn)為該SV是1,2,3 shared SV溶诞。
表格中的其他lineage分配類似鸯檬。
注意：
1）后期在分配時只考慮SV有沒有重疊，不需要考慮他們的坐標(biāo)螺垢，也不需要合并他們的坐標(biāo)喧务。
2）后期判斷某一個SV在物種中有沒有，就看他有沒有和物種里的SV重疊枉圃，只要有一個重疊功茴，就認(rèn)為它在該物種中有
1, AssignDeletion缺失
第一步：得到每兩兩物種之間overlap的SV列表，
只要兩個SV重疊度大于50%則認(rèn)為是同一個SV（這里采用的標(biāo)準(zhǔn)是參照reference【1】里的標(biāo)準(zhǔn)孽亲，也可以設(shè)置其他的標(biāo)準(zhǔn)）
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的做法進(jìn)行分配。
注：convert5SpeciesSVToLineageSV-v2.pl是小編自寫腳本返劲。如有需要玲昧，請進(jìn)QQ群下載，群號：756263353篮绿， 文件分享在群文件里面
參考資料：
【1】2018-Science- High-resolution comparative analysis of great apegenomes
