updates to support checking BAM headers match FASTA headers

Author: Jon Palmer

README.md

@@ -120,7 +120,7 @@ funannotate sort -i genome.cleaned.fa -o genome.final.fa
 hisat2-build genome.final.fa genome
 
 #now align reads to reference, using 12 cpus
-hisat2 --max_intronlen 3000 -p 12 -x genome -1 forward_R1.fastq -2 reverse_R2.fastq \
+hisat2 --max-intronlen 3000 -p 12 -x genome -1 forward_R1.fastq -2 reverse_R2.fastq \
        | samtools view -bS - | samtools sort -o RNA_seq.bam - 
 ```
 3b) You can also run something like Trinity/PASA/TransDecoder with the RNA-seq reads, for instructions see [Trinity](https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running%20Trinity) and [PASA](http://pasapipeline.github.io).

bin/funannotate-predict.py

@@ -215,6 +215,11 @@ def __init__(self, prog):
     if not header_test:
         lib.log.error("Fasta headers on your input have more characters than the max (%i), reformat headers to continue." % args.header_length)
         sys.exit(1)
+    #if BAM file passed, check if headers are same as input
+    if args.rna_bam:
+        if not lib.BamHeaderTest(args.input, args.rna_bam):
+            lib.log.error("Fasta headers in BAM file do not match genome, exiting.")
+            sys.exit(1)
     #just copy the input fasta to the misc folder and move on.
     shutil.copyfile(args.input, genome_input)
     Genome = os.path.abspath(genome_input)
@@ -224,6 +229,11 @@ def __init__(self, prog):
         sys.exit(1)
     for x in [args.masked_genome, args.repeatmasker_gff3]:
         lib.checkinputs(x)
+    #if BAM file passed, check if headers are same as input
+    if args.rna_bam:
+        if not lib.BamHeaderTest(args.masked_genome, args.rna_bam):
+            lib.log.error("Fasta headers in BAM file do not match genome, exiting.")
+            sys.exit(1)
     #now copy over masked genome and repeatmasker
     shutil.copyfile(args.masked_genome, genome_input)
     shutil.copyfile(args.masked_genome, MaskGenome)
@@ -239,7 +249,6 @@ def __init__(self, prog):
 Converter2 = os.path.join(EVM, 'EvmUtils', 'misc', 'augustus_GTF_to_EVM_GFF3.pl')
 EVM2proteins = os.path.join(EVM, 'EvmUtils', 'gff3_file_to_proteins.pl')
 
-
 #repeatmasker, run if not passed from command line
 if not os.path.isfile(MaskGenome):
     if not args.repeatmodeler_lib:

lib/library.py

@@ -440,6 +440,28 @@ def checkFastaHeaders(input, limit):
     else:
         return True
 
+def BamHeaderTest(genome, mapping):
+    import pybam
+    #get list of fasta headers from genome
+    genome_headers = []
+    with open(genome, 'rU') as input:
+        for rec in SeqIO.parse(input, 'fasta'):
+            if rec.id not in genome_headers:
+                genome_headers.append(rec.id)
+    #get list of fasta headers from BAM
+    bam_headers = []
+    with open(mapping, 'rU') as bamin:
+        bam_file = pybam.bgunzip(bamin)
+        bam_headers = bam_file.chromosomes_from_header
+    #now compare lists, basically if BAM headers not in genome headers, then output bad names to logfile and return FALSE
+    genome_headers = set(genome_headers)
+    diffs = [x for x in bam_headers if x not in genome_headers]
+    if len(diffs) > 0:
+        log.debug("ERROR: These BAM headers not found in genome FASTA headers\n%s" % ','.join(diffs))
+        return False
+    else:
+        return True
+    
 def gb2allout(input, GFF, Proteins, Transcripts, DNA):
     #this will not output any UTRs for gene models, don't think this is a problem right now....
     with open(GFF, 'w') as gff: