From c2f2cb03ced093ce87a4b3a661ba8c1f8faee4fa Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 6 Nov 2014 09:12:53 -0500 Subject: [PATCH] to be moved to a separate project --- bwa-typeHLA.sh | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/bwa-typeHLA.sh b/bwa-typeHLA.sh index b91067b..1cc74c4 100755 --- a/bwa-typeHLA.sh +++ b/bwa-typeHLA.sh @@ -21,18 +21,19 @@ if [ $is_ctg -eq 0 ]; then $root/fermi2.pl unitig -t2 -l$len -p $pre.tmp $pre.fq > $pre.tmp.mak make -f $pre.tmp.mak EXE_FERMI2=$root/fermi2 EXE_ROPEBWT2=$root/ropebwt2 else - ln -sf $pre.fq $pre.tmp.mag.gz + rm -f $pre.tmp.mag.gz + ln -s $pre.fq $pre.tmp.mag.gz fi # get contigs overlapping HLA exons (ls $root/HLA-idx/*.fa | xargs -i $root/bwa mem -t2 -B1 -O1 -E1 {} $pre.tmp.mag.gz 2>/dev/null) | grep -v ^@ | sort -k3,3 -k4,4n | bgzip > $pre.tmp.sam.gz $root/tabix -Bpsam $pre.tmp.sam.gz $root/HLA-approx.anno | cut -f1 | sort | uniq > $pre.tmp.kept -gzip -dc $pre.tmp.sam.gz | perl -ane 'print "$F[0]\n" if /AS:i:(\d+).*XS:i:(\d+)/&&$1==$2' | cut -f1 | sort | uniq > $pre.tmp.dropped +gzip -dc $pre.tmp.sam.gz | perl -ane 'print if $F[5]!~/^\d+[SH]/ && $F[5]!~/[SH]$/ && /AS:i:(\d+).*XS:i:(\d+)/ && $1==$2' | cut -f1 | sort | uniq > $pre.tmp.dropped awk -v f=$pre.tmp.dropped 'BEGIN{while((getline0)l[$1]=1}!l[$1]' $pre.tmp.kept | $root/seqtk subseq $pre.tmp.mag.gz - | gzip -1 > $pre.tmp.fq.gz # map HLA exons to de novo contigs -$root/bwa index -p $pre.tmp $pre.tmp.fq.gz 2> /dev/null -$root/bwa mem -aD.1 -t2 $pre.tmp $root/hla.fa | egrep "^(@|$2)" > $pre.sam +$root/bwa index -p $pre.tmp $pre.tmp.fq.gz 2>/dev/null +$root/bwa mem -aD.1 -t2 $pre.tmp $root/hla.fa 2>/dev/null | egrep "^(@|$2)" > $pre.sam # type HLA $root/k8 $root/bwa-typeHLA.js $pre.sam > $pre.gt