Update HBVIF.v0.2.1.sh

[DomeJoyce]

Run BLAT instead of BLAST

cat new_list_to_continue | while read line; do
let COUNTER+=1
echo "$line, which is the integration event $COUNTER out of $numberOfBlastRemaining"

Run BLAT with equivalent parameters

blat ../../genome/HG_HBV $line $line.psl -t=dna -q=dna -stepSize=7 -tileSize=7 -minScore=0 -repMatch=1024 -minIdentity=70 -mask=lower -maxIntron=500000 -minMatch=2 -minRepDivergence=0.001

Convert BLAT PSL output to BLAST-like output

awk 'BEGIN {OFS="\t"} {if ($0 !~ /^#/ && $0 !~ /psLayout/) {print $14,$16,$17,$10,$11}}' $line.psl > $line.intermediate.blat.out
done

Process the intermediate BLAT output to extract HBV hits

for line in *.intermediate.blat.out; do
grep "NC_003977.2" $line | head -n 1 > $line.HBV
done

Further process the intermediate BLAT output to extract human hits based on coordinates from BED file

for line in *.intermediate.blat.out; do
E=$(echo $line | cut -f1 -d ".")
C=$(awk -F '\t' -v E=$E '{split($4,a,";");if(a[1]==E) print $1}' ../../SA_primary_secondary_no_HBV.bed)
S=$(awk -F '\t' -v E=$E '{split($4,a,";");if(a[1]==E) print $2}' ../../SA_primary_secondary_no_HBV.bed)
F=$(awk -F '\t' -v E=$E '{split($4,a,";");if(a[1]==E) print $3}' ../../SA_primary_secondary_no_HBV.bed)
awk -F '\t' -v C="$C" -v S="$S" -v F="$F" '{if($1==C && $2>=S-100 && $3<=F+100) print $0}' $line | head -n1 > $line.HUMAN
done

Final processing of the BLAT output

cat list_for_paste_final | grep "out." | paste - - | while read line; do
awk '{split(FILENAME,a,"."); print $0 > a[1]"."a[2]".final.blat.out"}' $line
done

Convert each final BLAT output to a single line

for l in *.final.blat.out; do
paste -s $l > $l".one_line"
done

Move intermediate and final results to the appropriate directory

for f in *.blat.out; do mv $f ./intermediate_results/; done
for f in *.blat.out.HBV; do mv $f ./intermediate_results/; done
for f in *.blat.out.HUMAN; do mv $f ./intermediate_results/; done