zzh
/
bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

added the pipeline script

This commit is contained in:
Heng Li 2014-11-03 15:16:14 -05:00
parent 40c4a6ffe1
commit 9dbf152e31
2 changed files with 27 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
bwa-typeHLA.js
View File

@ -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];
}

19
bwa-typeHLA.sh 100755
View File

@ -0,0 +1,19 @@
#!/bin/bash
if [ $# -lt 2 ]; then
echo "Usage: $0 <prefix> <gene>"
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}