BLAST long sequences and get BAM files - wdlingit/cop GitHub Wiki

Temp page.

NOTE: should mention the limit on sequence header

Example command 1

$ ml BLAST+/2.13.0-gompi-2022a

$ makeblastdb -dbtype nucl -in NIES2134.fa

(blast TV2023 against NIES2134)
$ blastn -task blastn -query TV2023.fa -db NIES2134.fa > TV2023-NIES2134.comp.blastn.out

$ ml RackJ

$ preindex.pl NIES2134.fa

$ cat TV2023-NIES2134.comp.blastn.out | blast2psl.pl -mode blastn | psl2sam.pl -md NIES2134.fa /dev/stdin | samGetSEQ.pl -hardClip 1 /dev/stdin TV2023.fa | perl -ne 'if(/^\@/){ print }else{ @t=split; $headClip=0; $tailClip=0; if($t[5]=~/^(\d+)H(.+)/){ $headClip=$1; } if($t[5]=~/(.+?)(\d+)H$/){ $tailClip=$2; } if($t[1]&16){ $t[0].="_$.\_$tailClip" }else{ $t[0].="_$.\_$headClip" } print join("\t",@t)."\n" }' | java misc.AlignmentFilter2 -M SAM /dev/stdin -O /dev/stdout -quiet true -filter alnLen 150 -filter mID 0.85 | samtools view -Sbo /dev/stdout -T NIES2134.fa /dev/stdin | samtools calmd -bS /dev/stdin NIES2134.fa 2>/dev/null > TV2023-NIES2134.comp.blastn.bam

Example command 2

$ ls *.fa | perl -ne 'chomp; $cmd="formatdb -i $_ -p F"; print "$cmd\n"; system $cmd'

$ ls *.fa | perl -ne 'chomp; /(.+?)\./; $cmd="blastall -p blastn -i $_ -d $_ > $1.self.blastn.out"; print "$cmd\n"; system $cmd'

$ cat B1_cyano05_02.self.blastn.out | blast2psl.pl -mode blastn | psl2sam.pl -md B1_cyano05_02.fa /dev/stdin | perl -ne 'chomp; @t=split; if($t[5]=~/^(\d+)H/){ $clip=$1 }else{ $clip=0 } if( (($t[1]&16)==0) && ($t[0] eq $t[2]) && ($t[3]==($clip+1)) ){}else{ print "$_\n" }' | samGetSEQ.pl -hardClip 1 /dev/stdin B1_cyano05_02.fa | perl -ne 'if(/^\@/){ print }else{ @t=split; $headClip=0; $tailClip=0; if($t[5]=~/^(\d+)H(.+)/){ $headClip=$1; } if($t[5]=~/(.+?)(\d+)H$/){ $tailClip=$2; } if($t[1]&16){ $t[0].="_$.\_$tailClip" }else{ $t[0].="_$.\_$headClip" } print join("\t",@t)."\n" }' | java -classpath /home/wdlin/Tools/rackJ/rackj.jar:/home/bictools/Tools/sam-1.89.jar misc.AlignmentFilter2 -M SAM /dev/stdin -O /dev/stdout -quiet true -filter alnLen 150 -filter mID 0.85 | samtools view -Sbo /dev/stdout -T B1_cyano05_02.fa /dev/stdin | samtools calmd -bS /dev/stdin B1_cyano05_02.fa 2>/dev/null > B1_cyano05_02.self.blastn.bam