diff --git a/bwa-typeHLA.js b/bwa-typeHLA.js index 5561952..ad47387 100644 --- a/bwa-typeHLA.js +++ b/bwa-typeHLA.js @@ -195,14 +195,14 @@ for (var e = 0; e < exons.length; ++e) { } // drop mismapped contigs var dropped = []; - for (var c = 0; c < ca.length; ++c) { - var min = 0x7fffffff, cc = ca[c]; - for (var g = 0; g < ga.length; ++g) { - var gg = ga[g]; - min = min < sc[gg][cc]? min : sc[gg][cc]; + for (var k = 0; k < ca.length; ++k) { + var min = 0x7fffffff, c = ca[k]; + for (var i = 0; i < ga.length; ++i) { + var g = ga[i]; + min = min < sc[g][c]? min : sc[g][c]; } - dropped[cc] = min > drop_thres? true : false; - if (dropped[cc]) warn("Dropped contig " +clist[cc]+ " due to high divergence to all genes (minNM=" +min+ ")"); + dropped[c] = min > drop_thres? true : false; + if (dropped[c]) warn("Dropped contig " +clist[c]+ " due to high divergence to all genes (minNM=" +min+ ")"); } // fill the pair array var min_nm = 0xffff; @@ -211,7 +211,7 @@ for (var e = 0; e < exons.length; ++e) { for (var j = i; j < ga.length; ++j) { var gj = ga[j], g2 = sc[gj], m = 0; for (var k = 0; k < ca.length; ++k) { - c = ca[k]; + var c = ca[k]; if (!dropped[c]) m += g1[c] < g2[c]? g1[c] : g2[c]; } diff --git a/bwa-typeHLA.sh b/bwa-typeHLA.sh new file mode 100755 index 0000000..cc4bd12 --- /dev/null +++ b/bwa-typeHLA.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +if [ $# -lt 2 ]; then + echo "Usage: $0 " + exit 1 +fi + +root=`dirname $0` +pre=$1.$2 +len=`$root/seqtk comp $pre.fq | awk '{++x;y+=$2}END{printf("%.0f\n", y/x)}'` + +$root/fermi2.pl unitig -t2 -l$len -p $pre.utg $pre.fq > $pre.utg.mak +make -f $pre.utg.mak +$root/bwa mem -B1 -O1 -E1 $root/HLA-ALT.fa $pre.utg.mag.gz | grep -v ^@ | sort -k3,3 -k4,4n | bgzip > $pre.utg.alt.gz +$root/tabix -Bpsam $pre.utg.alt.gz $root/HLA-approx.anno | cut -f1 | sort | uniq | $root/seqtk subseq $pre.utg.mag.gz - | gzip -1 > $pre.utg.fq.gz +rm $pre.utg.{flt,ec,raw,pre}.* $pre.utg.mag* $pre.utg.mak $pre.utg.alt.gz +$root/bwa index -p $pre.utg $pre.utg.fq.gz +$root/bwa mem -aD.1 -t2 $pre.utg $root/hla.fa | egrep "^(@|$2)" > $pre.utg.sam +rm $pre.utg.{ann,amb,bwt,sa,pac}