ipyrad/assemble/refmap.py
def bam_region_to_fasta(data, sample, proc1, chrom, region_start, region_end):
    """ 
    Take the chromosome position, and start and end bases and return sequences
    of all reads that overlap these sites. This is the command we're building:

    samtools view -b 1A_sorted.bam 1:116202035-116202060 | \
             samtools bam2fq <options> -

            -b      : output bam format
            -0      : For SE, output all reads to this file
            -1/-2   : For PE, output first and second reads to different files
            -       : Tell samtools to read in from the pipe

    Write out the sam output and parse it to return as fasta for clust.gz file. 
    We also grab the reference sequence with a @REF header to aid in alignment
    for single-end data. This will be removed post-alignment. 
    """

    ## output bam file handle for storing genome regions
    bamf = sample.files.mapped_reads
    if not os.path.exists(bamf):
        raise IPyradWarningExit("  file not found - %s", bamf)

    # chrom = re.escape(repr(chrom))[1:-1].replace('\\\\', '\\')
    #LOGGER.info("before: %s", chrom)
    chrom.replace("|", r"\|")
    chrom.replace("(", r"\(")
    chrom.replace(")", r"\)")
    #LOGGER.info("after: %s", chrom)

    ## What we want to do is have the num-chrom dict as an arg, then build this
    ## string as three ints [chrom-int, pos-start, pos-end]
    #cint = cdict[chrom]
    #cpstring = "__{}_{}_{}__".format(cint, int(region_start)+1, region_end)

    ## a string argument as input to commands, indexed at either 0 or 1, 
    ## and with pipe characters removed from chromo names
    ## rstring_id1 is for fetching the reference sequence bcz faidx is
    ## 1 indexed
    rstring_id1 = "{}:{}-{}"\
        .format(chrom, str(int(region_start)+1), region_end)

    ## rstring_id0 is just for printing out the reference CHROM/POS
    ## in the read name
    rstring_id0 = "{}:{}-{}"\
        .format(chrom, region_start, region_end)

    ## If SE then we enforce the minimum overlap distance to avoid the 
    ## staircase syndrome of multiple reads overlapping just a little.
    overlap_buffer = 0
    if not "pair" in data.paramsdict["datatype"]:
        overlap_buffer = data._hackersonly["min_SE_refmap_overlap"]

    ## rstring_id0_buffered is for samtools view. We have to play patty
    ## cake here with the two rstring_id0s because we want `view` to 
    ## enforce the buffer for SE, but we want the reference sequence
    ## start and end positions to print correctly for downstream.
    rstring_id0_buffered = "{}:{}-{}"\
        .format(chrom, int(region_start) + overlap_buffer,\
                int(region_end) - overlap_buffer)

    ## The "samtools faidx" command will grab this region from reference 
    ## which we'll paste in at the top of each stack to aid alignment.
    cmd1 = [ipyrad.bins.samtools, "faidx", 
            data.paramsdict["reference_sequence"], 
            rstring_id1, " ; echo __done__"]

    ## Call the command, I found that it doesn't work with shell=False if 
    ## the refstring is 'MT':100-200', but it works if it is MT:100-200. 
    LOGGER.info("Grabbing bam_region_to_fasta:\n {}".format(cmd1))
    #proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE)
    #ref = proc1.communicate()[0]
    #if proc1.returncode:
    #    raise IPyradWarningExit("  error in %s: %s", cmd1, ref)

    ## push the samtools faidx command to our subprocess, then accumulate
    ## the results from stdout
    print(" ".join(cmd1), file=proc1.stdin)
    ref = ""
    for line in iter(proc1.stdout.readline, "//\n"):
        if "__done__" in line:
            break
        ref += line

    ## initialize the fasta list.
    fasta = []

    ## parse sam to fasta. Save ref location to name.
    ## Set size= an improbably large value so the REF sequence
    ## sorts to the top for muscle aligning.
    try:
        name, seq = ref.strip().split("\n", 1)
        seq = "".join(seq.split("\n"))
        fasta = ["{}_REF;size={};+\n{}".format(name, 1000000, seq)]
    except ValueError as inst:
        LOGGER.error("ref failed to parse - {}".format(ref))
        LOGGER.error(" ".join(cmd1))

    ## if PE then you have to merge the reads here
    if "pair" in data.paramsdict["datatype"]:
        ## PE R1 can either be on the forward or the reverse strand.
        ## Samtools view always outputs reads with respect to the
        ## forward strand. This means that reads with R1 on reverse
        ## end up with the R1 and R2 reference sequences swapped
        ## in the clust.gz file. There is a way to fix it but it's
        ## very annoying and i'm not sure if it's worth it...
        ## Drop the reference sequence for now...
        ##
        ## If you ever fix this be sure to remove the reference sequence
        ## from each cluster post alignment in cluster_within/align_and_parse()
        fasta = []

        ## Create temporary files for R1, R2 and merged, which we will pass to
        ## the function merge_pairs() which calls vsearch to test merging.
        ##
        ## If you are on linux then creating the temp files in /dev/shm
        ## should improve performance
        if os.path.exists("/dev/shm"):
            prefix = os.path.join("/dev/shm",
                            "{}-{}".format(sample.name, rstring_id0))
        else:
            prefix = os.path.join(data.dirs.refmapping, 
                            "{}-{}".format(sample.name, rstring_id0))
        read1 = "{}-R1".format(prefix)
        read2 = "{}-R2".format(prefix)
        merged = "{}-merged".format(prefix)

        ## Grab all the reads that map to this genomic location and dump
        ## fastq to R1 and R2 files.
        ## `-v 45` sets the default qscore to something high
        cmd1 = " ".join([ipyrad.bins.samtools, "view", "-b", bamf, rstring_id0])
        cmd2 = " ".join([ipyrad.bins.samtools, "bam2fq", "-v", "45", "-1", read1, "-2", read2, "-", "; echo __done__"])
        cmd = " | ".join([cmd1, cmd2])

        print(cmd, file=proc1.stdin)
        for line in iter(proc1.stdout.readline, "//\n"):
            if "__done__" in line:
                break

        ## run commands, pipe 1 -> 2, then cleanup
        ## proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE)
        ## proc2 = sps.Popen(cmd2, stderr=sps.STDOUT, stdout=sps.PIPE, stdin=proc1.stdout)
        ## res = proc2.communicate()[0]
        ## if proc2.returncode:
        ##     raise IPyradWarningExit("error {}: {}".format(cmd2, res))
        ## proc1.stdout.close()

        ## merge the pairs. 0 means don't revcomp bcz samtools already
        ## did it for us. 1 means "actually merge".
        try:
            ## return number of merged reads, writes merged data to 'merged'
            ## we don't yet do anything with the returned number of merged 
            _ = merge_pairs(data, [(read1, read2)], merged, 0, 1)

            with open(merged, 'r') as infile:
                quatro = itertools.izip(*[iter(infile)]*4)
                while 1:
                    try:
                        bits = quatro.next()
                    except StopIteration:
                        break
                    ## TODO: figure out a real way to get orientation for PE
                    orient = "+"
                    fullfast = ">{a};{b};{c};{d}\n{e}".format(
                        a=bits[0].split(";")[0],
                        b=rstring_id1,
                        c=bits[0].split(";")[1], 
                        d=orient,
                        e=bits[1].strip())
                    #,e=bits[9])
                    fasta.append(fullfast)

                ## TODO: If you ever figure out a good way to get the reference
                ## sequence included w/ PE then this commented call is useful
                ## for trimming the reference sequence to be the right length.
                ## If doing PE and R1/R2 don't overlap then the reference sequence
                ## will be quite long and will cause indel hell during the 
                ## alignment stage. Here trim the reference sequence to the length
                ## of the merged reads.
                ## This is commented out because we aren't currently including the
                ## ref seq for PE alignment.
                #fasta = trim_reference_sequence(fasta)

        except (OSError, ValueError, IPyradError) as inst:
            ## ValueError raised inside merge_pairs() if it can't open one
            ## or both of the files. Write this out, but ignore for now.
            ## Failed merging, probably unequal number of reads in R1 and R2
            ## IPyradError raised if merge_pairs can't read either R1 or R2
            ## file.
            ## Skip this locus?
            LOGGER.debug("Failed to merge reads, continuing; %s", inst)
            LOGGER.error("cmd - {}".format(cmd))
            return ""
        finally:
            ## Only clean up the files if they exist otherwise it'll raise.
            if os.path.exists(merged):
                os.remove(merged)
            if os.path.exists(read1):
                os.remove(read1)
            if os.path.exists(read2):
                os.remove(read2)
       
    else:
        try:
            ## SE if faster than PE cuz it skips writing intermedidate files
            ## rstring_id0_buffered is going to enforce the required 
            ## min_SE_refmap_overlap on either end of this region.
            cmd2 = [ipyrad.bins.samtools, "view", bamf, rstring_id0_buffered]
            proc2 = sps.Popen(cmd2, stderr=sps.STDOUT, stdout=sps.PIPE)
    
            ## run and check outputs
            res = proc2.communicate()[0]
            if proc2.returncode:
                raise IPyradWarningExit("{} {}".format(cmd2, res))

            ## if the region string is malformated you'll get back a warning
            ## from samtools
            if "[main_samview]" in res:
                raise IPyradError("Bad reference region {}".format(rstring_id0_buffered))

            ## do not join seqs that
            for line in res.strip().split("\n"):
                bits = line.split("\t")

                ## Read in the 2nd field (FLAGS), convert to binary
                ## and test if the 7th bit is set which indicates revcomp
                orient = "+"
                if int('{0:012b}'.format(int(bits[1]))[7]):
                    orient = "-"
                    ## Don't actually revcomp the sequence because samtools
                    ## writes out reference sequence on the forward strand
                    ## as well as reverse strand hits from the bam file.
                    #bits[9] = revcomp(bits[9])

                ## Rip insert the mapping position between the seq label and
                ## the vsearch derep size.
                fullfast = ">{a};{b};{c};{d}\n{e}".format(
                    a=bits[0].split(";")[0],
                    b=rstring_id0,
                    c=bits[0].split(";")[1],
                    d=orient,
                    e=bits[9])
                fasta.append(fullfast)
        except IPyradError as inst:
            ## If the mapped fragment is too short then the you'll get
            ## regions that look like this: scaffold262:299039-299036
            ## Just carry on, it's not a big deal.
            LOGGER.debug("Got a bad region string: {}".format(inst))
            return ""
        except (OSError, ValueError, Exception) as inst:
            ## Once in a blue moon something fsck and it breaks the
            ## assembly. No reason to give up if .001% of reads fail
            ## so just skip this locus.
            LOGGER.error("Failed get reads at a locus, continuing; %s", inst)
            LOGGER.error("cmd - {}".format(cmd2))
            return ""

    return "\n".join(fasta)
def bam_region_to_fasta(data, sample, proc1, chrom, region_start, region_end):
    """ 
    Take the chromosome position, and start and end bases and return sequences
    of all reads that overlap these sites. This is the command we're building:

    samtools view -b 1A_sorted.bam 1:116202035-116202060 | \
             samtools bam2fq <options> -

            -b      : output bam format
            -0      : For SE, output all reads to this file
            -1/-2   : For PE, output first and second reads to different files
            -       : Tell samtools to read in from the pipe

    Write out the sam output and parse it to return as fasta for clust.gz file. 
    We also grab the reference sequence with a @REF header to aid in alignment
    for single-end data. This will be removed post-alignment. 
    """

    ## output bam file handle for storing genome regions
    bamf = sample.files.mapped_reads
    if not os.path.exists(bamf):
        raise IPyradWarningExit("  file not found - %s", bamf)

    # chrom = re.escape(repr(chrom))[1:-1].replace('\\\\', '\\')
    #LOGGER.info("before: %s", chrom)
    chrom.replace("|", r"\|")
    chrom.replace("(", r"\(")
    chrom.replace(")", r"\)")
    #LOGGER.info("after: %s", chrom)

    ## What we want to do is have the num-chrom dict as an arg, then build this
    ## string as three ints [chrom-int, pos-start, pos-end]
    #cint = cdict[chrom]
    #cpstring = "__{}_{}_{}__".format(cint, int(region_start)+1, region_end)

    ## a string argument as input to commands, indexed at either 0 or 1, 
    ## and with pipe characters removed from chromo names
    ## rstring_id1 is for fetching the reference sequence bcz faidx is
    ## 1 indexed
    rstring_id1 = "{}:{}-{}"\
        .format(chrom, str(int(region_start)+1), region_end)

    ## rstring_id0 is just for printing out the reference CHROM/POS
    ## in the read name
    rstring_id0 = "{}:{}-{}"\
        .format(chrom, region_start, region_end)

    ## If SE then we enforce the minimum overlap distance to avoid the 
    ## staircase syndrome of multiple reads overlapping just a little.
    overlap_buffer = 0
    if not "pair" in data.paramsdict["datatype"]:
        overlap_buffer = data._hackersonly["min_SE_refmap_overlap"]

    ## rstring_id0_buffered is for samtools view. We have to play patty
    ## cake here with the two rstring_id0s because we want `view` to 
    ## enforce the buffer for SE, but we want the reference sequence
    ## start and end positions to print correctly for downstream.
    rstring_id0_buffered = "{}:{}-{}"\
        .format(chrom, int(region_start) + overlap_buffer,\
                int(region_end) - overlap_buffer)

    ## The "samtools faidx" command will grab this region from reference 
    ## which we'll paste in at the top of each stack to aid alignment.
    cmd1 = [ipyrad.bins.samtools, "faidx", 
            data.paramsdict["reference_sequence"], 
            rstring_id1, " ; echo __done__"]

    ## Call the command, I found that it doesn't work with shell=False if 
    ## the refstring is 'MT':100-200', but it works if it is MT:100-200. 
    LOGGER.info("Grabbing bam_region_to_fasta:\n {}".format(cmd1))
    #proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE)
    #ref = proc1.communicate()[0]
    #if proc1.returncode:
    #    raise IPyradWarningExit("  error in %s: %s", cmd1, ref)

    ## push the samtools faidx command to our subprocess, then accumulate
    ## the results from stdout
    print(" ".join(cmd1), file=proc1.stdin)
    ref = ""
    for line in iter(proc1.stdout.readline, "//\n"):
        if "__done__" in line:
            break
        ref += line

    ## initialize the fasta list.
    fasta = []

    ## parse sam to fasta. Save ref location to name.
    ## Set size= an improbably large value so the REF sequence
    ## sorts to the top for muscle aligning.
    try:
        name, seq = ref.strip().split("\n", 1)
        seq = "".join(seq.split("\n"))
        fasta = ["{}_REF;size={};+\n{}".format(name, 1000000, seq)]
    except ValueError as inst:
        LOGGER.error("ref failed to parse - {}".format(ref))
        LOGGER.error(" ".join(cmd1))

    ## if PE then you have to merge the reads here
    if "pair" in data.paramsdict["datatype"]:
        ## PE R1 can either be on the forward or the reverse strand.
        ## Samtools view always outputs reads with respect to the
        ## forward strand. This means that reads with R1 on reverse
        ## end up with the R1 and R2 reference sequences swapped
        ## in the clust.gz file. There is a way to fix it but it's
        ## very annoying and i'm not sure if it's worth it...
        ## Drop the reference sequence for now...
        ##
        ## If you ever fix this be sure to remove the reference sequence
        ## from each cluster post alignment in cluster_within/align_and_parse()
        fasta = []

        ## Create temporary files for R1, R2 and merged, which we will pass to
        ## the function merge_pairs() which calls vsearch to test merging.
        ##
        ## If you are on linux then creating the temp files in /dev/shm
        ## should improve performance
        if os.path.exists("/dev/shm"):
            prefix = os.path.join("/dev/shm",
                            "{}-{}".format(sample.name, rstring_id0))
        else:
            prefix = os.path.join(data.dirs.refmapping, 
                            "{}-{}".format(sample.name, rstring_id0))
        read1 = "{}-R1".format(prefix)
        read2 = "{}-R2".format(prefix)
        merged = "{}-merged".format(prefix)

        ## Grab all the reads that map to this genomic location and dump
        ## fastq to R1 and R2 files.
        ## `-v 45` sets the default qscore to something high
        cmd1 = " ".join([ipyrad.bins.samtools, "view", "-b", bamf, rstring_id0])
        cmd2 = " ".join([ipyrad.bins.samtools, "bam2fq", "-v", "45", "-1", read1, "-2", read2, "-", "; echo __done__"])
        cmd = " | ".join([cmd1, cmd2])

        print(cmd, file=proc1.stdin)
        for line in iter(proc1.stdout.readline, "//\n"):
            if "__done__" in line:
                break

        ## run commands, pipe 1 -> 2, then cleanup
        ## proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE)
        ## proc2 = sps.Popen(cmd2, stderr=sps.STDOUT, stdout=sps.PIPE, stdin=proc1.stdout)
        ## res = proc2.communicate()[0]
        ## if proc2.returncode:
        ##     raise IPyradWarningExit("error {}: {}".format(cmd2, res))
        ## proc1.stdout.close()

        ## merge the pairs. 0 means don't revcomp bcz samtools already
        ## did it for us. 1 means "actually merge".
        try:
            ## return number of merged reads, writes merged data to 'merged'
            ## we don't yet do anything with the returned number of merged 
            _ = merge_pairs(data, [(read1, read2)], merged, 0, 1)

            with open(merged, 'r') as infile:
                quatro = itertools.izip(*[iter(infile)]*4)
                while 1:
                    try:
                        bits = quatro.next()
                    except StopIteration:
                        break
                    ## TODO: figure out a real way to get orientation for PE
                    orient = "+"
                    fullfast = ">{a};{b};{c};{d}\n{e}".format(
                        a=bits[0].split(";")[0],
                        b=rstring_id1,
                        c=bits[0].split(";")[1], 
                        d=orient,
                        e=bits[1].strip())
                    #,e=bits[9])
                    fasta.append(fullfast)

                ## TODO: If you ever figure out a good way to get the reference
                ## sequence included w/ PE then this commented call is useful
                ## for trimming the reference sequence to be the right length.
                ## If doing PE and R1/R2 don't overlap then the reference sequence
                ## will be quite long and will cause indel hell during the 
                ## alignment stage. Here trim the reference sequence to the length
                ## of the merged reads.
                ## This is commented out because we aren't currently including the
                ## ref seq for PE alignment.
                #fasta = trim_reference_sequence(fasta)

        except (OSError, ValueError, IPyradError) as inst:
            ## ValueError raised inside merge_pairs() if it can't open one
            ## or both of the files. Write this out, but ignore for now.
            ## Failed merging, probably unequal number of reads in R1 and R2
            ## IPyradError raised if merge_pairs can't read either R1 or R2
            ## file.
            ## Skip this locus?
            LOGGER.debug("Failed to merge reads, continuing; %s", inst)
            LOGGER.error("cmd - {}".format(cmd))
            return ""
        finally:
            ## Only clean up the files if they exist otherwise it'll raise.
            if os.path.exists(merged):
                os.remove(merged)
            if os.path.exists(read1):
                os.remove(read1)
            if os.path.exists(read2):
                os.remove(read2)
       
    else:
        try:
            ## SE if faster than PE cuz it skips writing intermedidate files
            ## rstring_id0_buffered is going to enforce the required 
            ## min_SE_refmap_overlap on either end of this region.
            cmd2 = [ipyrad.bins.samtools, "view", bamf, rstring_id0_buffered]
            proc2 = sps.Popen(cmd2, stderr=sps.STDOUT, stdout=sps.PIPE)
    
            ## run and check outputs
            res = proc2.communicate()[0]
            if proc2.returncode:
                raise IPyradWarningExit("{} {}".format(cmd2, res))

            ## if the region string is malformated you'll get back a warning
            ## from samtools
            if "[main_samview]" in res:
                raise IPyradError("Bad reference region {}".format(rstring_id0_buffered))

            ## do not join seqs that
            for line in res.strip().split("\n"):
                bits = line.split("\t")

                ## Read in the 2nd field (FLAGS), convert to binary
                ## and test if the 7th bit is set which indicates revcomp
                orient = "+"
                if int('{0:012b}'.format(int(bits[1]))[7]):
                    orient = "-"
                    ## Don't actually revcomp the sequence because samtools
                    ## writes out reference sequence on the forward strand
                    ## as well as reverse strand hits from the bam file.
                    #bits[9] = revcomp(bits[9])

                ## Rip insert the mapping position between the seq label and
                ## the vsearch derep size.
                fullfast = ">{a};{b};{c};{d}\n{e}".format(
                    a=bits[0].split(";")[0],
                    b=rstring_id0,
                    c=bits[0].split(";")[1],
                    d=orient,
                    e=bits[9])
                fasta.append(fullfast)
        except IPyradError as inst:
            ## If the mapped fragment is too short then the you'll get
            ## regions that look like this: scaffold262:299039-299036
            ## Just carry on, it's not a big deal.
            LOGGER.debug("Got a bad region string: {}".format(inst))
            return ""
        except (OSError, ValueError, Exception) as inst:
            ## Once in a blue moon something fsck and it breaks the
            ## assembly. No reason to give up if .001% of reads fail
            ## so just skip this locus.
            LOGGER.error("Failed get reads at a locus, continuing; %s", inst)
            LOGGER.error("cmd - {}".format(cmd2))
            return ""

    return "\n".join(fasta)
def bam_region_to_fasta ( data , sample , proc1 , chrom , region_start , region_end ) : ## output bam file handle for storing genome regions bamf = sample . files . mapped_reads if not os . path . exists ( bamf ) : raise IPyradWarningExit ( "  file not found - %s" , bamf ) # chrom = re.escape(repr(chrom))[1:-1].replace('\\\\', '\\') #LOGGER.info("before: %s", chrom) chrom . replace ( "|" , r"\|" ) chrom . replace ( "(" , r"\(" ) chrom . replace ( ")" , r"\)" ) #LOGGER.info("after: %s", chrom) ## What we want to do is have the num-chrom dict as an arg, then build this ## string as three ints [chrom-int, pos-start, pos-end] #cint = cdict[chrom] #cpstring = "__{}_{}_{}__".format(cint, int(region_start)+1, region_end) ## a string argument as input to commands, indexed at either 0 or 1,  ## and with pipe characters removed from chromo names ## rstring_id1 is for fetching the reference sequence bcz faidx is ## 1 indexed rstring_id1 = "{}:{}-{}" . format ( chrom , str ( int ( region_start ) + 1 ) , region_end ) ## rstring_id0 is just for printing out the reference CHROM/POS ## in the read name rstring_id0 = "{}:{}-{}" . format ( chrom , region_start , region_end ) ## If SE then we enforce the minimum overlap distance to avoid the  ## staircase syndrome of multiple reads overlapping just a little. overlap_buffer = 0 if not "pair" in data . paramsdict [ "datatype" ] : overlap_buffer = data . _hackersonly [ "min_SE_refmap_overlap" ] ## rstring_id0_buffered is for samtools view. We have to play patty ## cake here with the two rstring_id0s because we want `view` to  ## enforce the buffer for SE, but we want the reference sequence ## start and end positions to print correctly for downstream. rstring_id0_buffered = "{}:{}-{}" . format ( chrom , int ( region_start ) + overlap_buffer , int ( region_end ) - overlap_buffer ) ## The "samtools faidx" command will grab this region from reference  ## which we'll paste in at the top of each stack to aid alignment. cmd1 = [ ipyrad . bins . samtools , "faidx" , data . paramsdict [ "reference_sequence" ] , rstring_id1 , " ; echo __done__" ] ## Call the command, I found that it doesn't work with shell=False if  ## the refstring is 'MT':100-200', but it works if it is MT:100-200.  LOGGER . info ( "Grabbing bam_region_to_fasta:\n {}" . format ( cmd1 ) ) #proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE) #ref = proc1.communicate()[0] #if proc1.returncode: #    raise IPyradWarningExit("  error in %s: %s", cmd1, ref) ## push the samtools faidx command to our subprocess, then accumulate ## the results from stdout print ( " " . join ( cmd1 ) , file = proc1 . stdin ) ref = "" for line in iter ( proc1 . stdout . readline , "//\n" ) : if "__done__" in line : break ref += line ## initialize the fasta list. fasta = [ ] ## parse sam to fasta. Save ref location to name. ## Set size= an improbably large value so the REF sequence ## sorts to the top for muscle aligning. try : name , seq = ref . strip ( ) . split ( "\n" , 1 ) seq = "" . join ( seq . split ( "\n" ) ) fasta = [ "{}_REF;size={};+\n{}" . format ( name , 1000000 , seq ) ] except ValueError as inst : LOGGER . error ( "ref failed to parse - {}" . format ( ref ) ) LOGGER . error ( " " . join ( cmd1 ) ) ## if PE then you have to merge the reads here if "pair" in data . paramsdict [ "datatype" ] : ## PE R1 can either be on the forward or the reverse strand. ## Samtools view always outputs reads with respect to the ## forward strand. This means that reads with R1 on reverse ## end up with the R1 and R2 reference sequences swapped ## in the clust.gz file. There is a way to fix it but it's ## very annoying and i'm not sure if it's worth it... ## Drop the reference sequence for now... ## ## If you ever fix this be sure to remove the reference sequence ## from each cluster post alignment in cluster_within/align_and_parse() fasta = [ ] ## Create temporary files for R1, R2 and merged, which we will pass to ## the function merge_pairs() which calls vsearch to test merging. ## ## If you are on linux then creating the temp files in /dev/shm ## should improve performance if os . path . exists ( "/dev/shm" ) : prefix = os . path . join ( "/dev/shm" , "{}-{}" . format ( sample . name , rstring_id0 ) ) else : prefix = os . path . join ( data . dirs . refmapping , "{}-{}" . format ( sample . name , rstring_id0 ) ) read1 = "{}-R1" . format ( prefix ) read2 = "{}-R2" . format ( prefix ) merged = "{}-merged" . format ( prefix ) ## Grab all the reads that map to this genomic location and dump ## fastq to R1 and R2 files. ## `-v 45` sets the default qscore to something high cmd1 = " " . join ( [ ipyrad . bins . samtools , "view" , "-b" , bamf , rstring_id0 ] ) cmd2 = " " . join ( [ ipyrad . bins . samtools , "bam2fq" , "-v" , "45" , "-1" , read1 , "-2" , read2 , "-" , "; echo __done__" ] ) cmd = " | " . join ( [ cmd1 , cmd2 ] ) print ( cmd , file = proc1 . stdin ) for line in iter ( proc1 . stdout . readline , "//\n" ) : if "__done__" in line : break ## run commands, pipe 1 -> 2, then cleanup ## proc1 = sps.Popen(cmd1, stderr=sps.STDOUT, stdout=sps.PIPE) ## proc2 = sps.Popen(cmd2, stderr=sps.STDOUT, stdout=sps.PIPE, stdin=proc1.stdout) ## res = proc2.communicate()[0] ## if proc2.returncode: ##     raise IPyradWarningExit("error {}: {}".format(cmd2, res)) ## proc1.stdout.close() ## merge the pairs. 0 means don't revcomp bcz samtools already ## did it for us. 1 means "actually merge". try : ## return number of merged reads, writes merged data to 'merged' ## we don't yet do anything with the returned number of merged  _ = merge_pairs ( data , [ ( read1 , read2 ) ] , merged , 0 , 1 ) with open ( merged , 'r' ) as infile : quatro = itertools . izip ( * [ iter ( infile ) ] * 4 ) while 1 : try : bits = quatro . next ( ) except StopIteration : break ## TODO: figure out a real way to get orientation for PE orient = "+" fullfast = ">{a};{b};{c};{d}\n{e}" . format ( a = bits [ 0 ] . split ( ";" ) [ 0 ] , b = rstring_id1 , c = bits [ 0 ] . split ( ";" ) [ 1 ] , d = orient , e = bits [ 1 ] . strip ( ) ) #,e=bits[9]) fasta . append ( fullfast ) ## TODO: If you ever figure out a good way to get the reference ## sequence included w/ PE then this commented call is useful ## for trimming the reference sequence to be the right length. ## If doing PE and R1/R2 don't overlap then the reference sequence ## will be quite long and will cause indel hell during the  ## alignment stage. Here trim the reference sequence to the length ## of the merged reads. ## This is commented out because we aren't currently including the ## ref seq for PE alignment. #fasta = trim_reference_sequence(fasta) except ( OSError , ValueError , IPyradError ) as inst : ## ValueError raised inside merge_pairs() if it can't open one ## or both of the files. Write this out, but ignore for now. ## Failed merging, probably unequal number of reads in R1 and R2 ## IPyradError raised if merge_pairs can't read either R1 or R2 ## file. ## Skip this locus? LOGGER . debug ( "Failed to merge reads, continuing; %s" , inst ) LOGGER . error ( "cmd - {}" . format ( cmd ) ) return "" finally : ## Only clean up the files if they exist otherwise it'll raise. if os . path . exists ( merged ) : os . remove ( merged ) if os . path . exists ( read1 ) : os . remove ( read1 ) if os . path . exists ( read2 ) : os . remove ( read2 ) else : try : ## SE if faster than PE cuz it skips writing intermedidate files ## rstring_id0_buffered is going to enforce the required  ## min_SE_refmap_overlap on either end of this region. cmd2 = [ ipyrad . bins . samtools , "view" , bamf , rstring_id0_buffered ] proc2 = sps . Popen ( cmd2 , stderr = sps . STDOUT , stdout = sps . PIPE ) ## run and check outputs res = proc2 . communicate ( ) [ 0 ] if proc2 . returncode : raise IPyradWarningExit ( "{} {}" . format ( cmd2 , res ) ) ## if the region string is malformated you'll get back a warning ## from samtools if "[main_samview]" in res : raise IPyradError ( "Bad reference region {}" . format ( rstring_id0_buffered ) ) ## do not join seqs that for line in res . strip ( ) . split ( "\n" ) : bits = line . split ( "\t" ) ## Read in the 2nd field (FLAGS), convert to binary ## and test if the 7th bit is set which indicates revcomp orient = "+" if int ( '{0:012b}' . format ( int ( bits [ 1 ] ) ) [ 7 ] ) : orient = "-" ## Don't actually revcomp the sequence because samtools ## writes out reference sequence on the forward strand ## as well as reverse strand hits from the bam file. #bits[9] = revcomp(bits[9]) ## Rip insert the mapping position between the seq label and ## the vsearch derep size. fullfast = ">{a};{b};{c};{d}\n{e}" . format ( a = bits [ 0 ] . split ( ";" ) [ 0 ] , b = rstring_id0 , c = bits [ 0 ] . split ( ";" ) [ 1 ] , d = orient , e = bits [ 9 ] ) fasta . append ( fullfast ) except IPyradError as inst : ## If the mapped fragment is too short then the you'll get ## regions that look like this: scaffold262:299039-299036 ## Just carry on, it's not a big deal. LOGGER . debug ( "Got a bad region string: {}" . format ( inst ) ) return "" except ( OSError , ValueError , Exception ) as inst : ## Once in a blue moon something fsck and it breaks the ## assembly. No reason to give up if .001% of reads fail ## so just skip this locus. LOGGER . error ( "Failed get reads at a locus, continuing; %s" , inst ) LOGGER . error ( "cmd - {}" . format ( cmd2 ) ) return "" return "\n" . join ( fasta )
Take the chromosome position, and start and end bases and return sequences
    of all reads that overlap these sites. This is the command we're building:

    samtools view -b 1A_sorted.bam 1:116202035-116202060 | \
             samtools bam2fq <options> -

            -b      : output bam format
            -0      : For SE, output all reads to this file
            -1/-2   : For PE, output first and second reads to different files
            -       : Tell samtools to read in from the pipe

    Write out the sam output and parse it to return as fasta for clust.gz file. 
    We also grab the reference sequence with a @REF header to aid in alignment
    for single-end data. This will be removed post-alignment.
Take the chromosome position and start and end bases and return sequences of all reads that overlap these sites . This is the command we re building :