<commit_before>#Ancient DNA pulldown - pulldown from bam into pseudo-haploid data in vcf format, or readcount output
#Python3... 

from __future__ import division, print_function
import argparse, gdc, pysam, random, pdb
from scipy.stats import binom

LOG10=2.302585092994046

################################################################################

def parse_options():
    """
    argparse
    """
    parser=argparse.ArgumentParser()
    parser.add_argument('-b', '--bamlist', type=str, default="", help=
                        "list of bam files. Two columns, with sample name in first column and bam file in second column")
    parser.add_argument('-s', '--snps', type=str, default="", help=
                        "File with snps to include, either in eigenstrat .snp or .vcf format, must be sorted")
    parser.add_argument('-o', '--output', type=str, default="vcf", help=
                        "Output type, either vcf or readcounts")
    parser.add_argument('-m', '--method', type=str, default="random", help=
                        "Output type, either random or majority")
    parser.add_argument('--seed', type=int, default=12345, help=
                        "Random number generator seed")
    parser.add_argument('-e', '--error', type=float, default=0.01, help=
                        "Error rate (for likelihood calculations)")
    parser.add_argument('-d', '--damage', type=float, default=0.05, help=
                        "Damage rate (for likelihood calculations)")
    parser.add_argument('-q', '--basequal', type=int, default=30, help=
                        "Minimum base quality")
    parser.add_argument('-i', '--mapqual', type=int, default=0, help=
                        "Minimum map quality")
    parser.add_argument('--gl', dest='gl', action='store_true')
    parser.set_defaults(likelihoods=False)


    args=parser.parse_args()
    if args.output not in ["vcf", "readcounts"]:
        raise Exception("output must be either vcf or readcounts")
    if args.method not in ["random", "majority"]:
        raise Exception("output must be either random or majority")

    return parser.parse_args()

################################################################################

class bamcollection:
    """
    Hold all the bamfiles and access the pileups
    """

    def __init__(self, bamlist, options):
        self.sample_names=[]
        self.bamfiles=[]
        self.min_base_qual=options.basequal
        self.min_map_qual=options.mapqual

        bams=open(bamlist)
        for line in bams:
            bits=line.split()
            self.sample_names.append(bits[0])
            self.bamfiles.append(pysam.AlignmentFile(bits[1],"rb"))
            
    def bases(self, chrom, pos):
        """
        For each bam, return the bases at each position. 
        """
        bases=[]
        for bf in self.bamfiles:
            base=[]
            for pupcol in bf.pileup( chrom, pos-1, pos):
                for pupread in pupcol.pileups: 
                    if pupcol.pos == pos-1 and pupread.query_position !=None:
                        bq=pupread.alignment.query_qualities[pupread.query_position]>self.min_base_qual
                        mq=pupread.alignment.mapq>self.min_map_qual
                        if bq and mq:
                            base.append(pupread.alignment.query_sequence[pupread.query_position])
            bases.append("".join(base))
        return(bases)

################################################################################

class snplist:
    """
    An iterator over a list of snps
    """
    
    def __init__(self, file):
        self.file=file
        if file[-4:]==".snp":
            self.is_snp=True
            self.is_vcf=False
            self.data=gdc.open2(file)
        elif file[-4:]==".vcf" or file[-7:]==".vcf.gz":
            self.is_snp=False<commit_msg>fix_merge<commit_after>#Ancient DNA pulldown - pulldown from bam into pseudo-haploid data in vcf format, or readcount output
#Python3... 

from __future__ import division, print_function
import argparse, gdc, pysam, random, pdb
from scipy.stats import binom

LOG10=2.302585092994046

################################################################################

def parse_options():
    """
    argparse
    """
    parser=argparse.ArgumentParser()
    parser.add_argument('-b', '--bamlist', type=str, default="", help=
                        "list of bam files. Two columns, with sample name in first column and bam file in second column")
    parser.add_argument('-s', '--snps', type=str, default="", help=
                        "File with snps to include, either in eigenstrat .snp or .vcf format, must be sorted")
    parser.add_argument('-o', '--output', type=str, default="vcf", help=
                        "Output type, either vcf or readcounts")
    parser.add_argument('-m', '--method', type=str, default="random", help=
                        "Output type, either random or majority")
    parser.add_argument('--seed', type=int, default=12345, help=
                        "Random number generator seed")
    parser.add_argument('-e', '--error', type=float, default=0.01, help=
                        "Error rate (for likelihood calculations)")
    parser.add_argument('-d', '--damage', type=float, default=0.05, help=
                        "Damage rate (for likelihood calculations)")
    parser.add_argument('-q', '--basequal', type=int, default=30, help=
                        "Minimum base quality")
    parser.add_argument('-i', '--mapqual', type=int, default=0, help=
                        "Minimum map quality")
    parser.add_argument('--gl', dest='gl', action='store_true')
    parser.set_defaults(likelihoods=False)


    args=parser.parse_args()
    if args.output not in ["vcf", "readcounts"]:
        raise Exception("output must be either vcf or readcounts")
    if args.method not in ["random", "majority"]:
        raise Exception("output must be either random or majority")

    return parser.parse_args()

################################################################################

class bamcollection:
    """
    Hold all the bamfiles and access the pileups
    """

    def __init__(self, bamlist, options):
        self.sample_names=[]
        self.bamfiles=[]
        self.min_base_qual=options.basequal
        self.min_map_qual=options.mapqual

        bams=open(bamlist)
        for line in bams:
            bits=line.split()
            self.sample_names.append(bits[0])
            self.bamfiles.append(pysam.AlignmentFile(bits[1],"rb"))
            
    def bases(self, chrom, pos):
        """
        For each bam, return the bases at each position. 
        """
        bases=[]
        for bf in self.bamfiles:
            base=[]

            for pupcol in bf.pileup( chrom, pos-1, pos):
                for pupread in pupcol.pileups: 
                    if pupcol.pos == pos-1 and pupread.query_position !=None:
                        bq=pupread.alignment.query_qualities[pupread.query_position]>=self.min_base_qual
                        mq=pupread.alignment.mapq>=self.min_map_qual
                        if bq and mq:
                            base.append(pupread.alignment.query_sequence[pupread.query_position])
            bases.append("".join(base))
        return(bases)

################################################################################

class snplist:
    """
    An iterator over a list of snps
    """
    
    def __init__(self, file):
        self.file=file
        if file[-4:]==".snp":
            self.is_snp=True
            self.is_vcf=False
            self.data=gdc.open2(file)
        elif file[-4:]==".vcf" or file[-7:]==".vcf.gz":
            self.is_snp=False
{'fraction_non_alphanumeric': 0.17804154302670624, 'fraction_numerical': 0.01057928009289124, 'mean_word_length': 2.391076115485564, 'pattern_counts': {'":': 4, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1243360', 'n_tokens_mistral': 2306, 'n_tokens_neox': 2187, 'n_words': 547}