diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e443ffb --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +/test_data +/build +/ATaRVa.egg-info +/ATARVA/__pycache__ +.DS_store \ No newline at end of file diff --git a/ATARVA/amplicon_read.py b/ATARVA/amplicon_read.py new file mode 100644 index 0000000..f282304 --- /dev/null +++ b/ATARVA/amplicon_read.py @@ -0,0 +1,24 @@ + + +def check_flank(cooper, read, locus_start, locus_end): + flank = cooper.args.flank + enough_soft_flank = False + left_flank = None; right_flank = None + + # read alignment covers the locus with enough soft clipped flanks on both sides + if read.ref_start < locus_start - flank and read.ref_end > locus_end + flank: + enough_soft_flank = True + left_flank = flank + right_flank = flank + + # for upstream soft_clip, checking the length compared to flank length + elif locus_start < read.ref_start: + mod_locus_start = read.ref_start + left_flank = min(flank, locus_start - read.ref_start) + + # for downstream soft_clip, checking the length compared to flank length + elif locus_end > read.ref_end: + mod_locus_end = read.ref_end + right_flank = min(flank, read.ref_end - locus_end) + + return enough_soft_flank, left_flank, right_flank \ No newline at end of file diff --git a/ATARVA/baseline.py b/ATARVA/baseline.py index 4772ae3..8b2dfe4 100644 --- a/ATARVA/baseline.py +++ b/ATARVA/baseline.py @@ -1,508 +1,525 @@ -from ATARVA.locus_utils import process_locus -from ATARVA.cstag_utils import parse_cstag -from ATARVA.cigar_utils import parse_cigar_tag -from ATARVA.operation_utils import update_homopolymer_coords -from ATARVA.genotype_utils import analyse_genotype -from ATARVA.vcf_writer import * - +import numpy as np import pysam -import sys -import logging - -def locus_processor(global_loci_keys, global_loci_ends, global_loci_variations, global_read_variations, global_snp_positions, prev_reads, sorted_global_snp_list, maxR, minR, ref, Chrom, global_loci_info, out, snpQ, snpC, snpD, snpR, phasingR, tbx, flank, sorted_global_ins_rpos_set, log_bool, logger, male): - - genotyped_loci = 0 - popped = global_loci_ends.pop(0) - locus_key = global_loci_keys.pop(0) - lstart = int(locus_key[locus_key.index(':')+1 : locus_key.index('-')]) - lend = int(locus_key[locus_key.index('-')+1:]) - near_by_loci = [] - for row in tbx.fetch(Chrom, lstart-flank, lend+flank): - row = row.split('\t') - near_by_loci.append( ( int(row[1]), int(row[2]) ) ) - - if locus_key in global_loci_variations: - - - if sorted_global_snp_list == []: - sorted_global_snp_list = sorted(list(global_snp_positions.keys())) - - - prev_reads, homozygous, ambiguous, homozygous_allele, reads_of_homozygous, hallele_counter, skip_point, max_limit = process_locus(locus_key, global_loci_variations, global_read_variations, global_snp_positions, prev_reads, sorted_global_snp_list, maxR, minR, global_loci_info, near_by_loci, sorted_global_ins_rpos_set, Chrom, lstart, lend, ref, log_bool, logger) - - if homozygous: - vcf_homozygous_writer(ref, Chrom, locus_key, global_loci_info, homozygous_allele, global_loci_variations, len(reads_of_homozygous), out, reads_of_homozygous) - genotyped_loci += 1 - elif ambiguous: - state, skip_point = analyse_genotype(Chrom, locus_key, global_loci_info, global_loci_variations, global_read_variations, global_snp_positions, hallele_counter, ref, out, sorted_global_snp_list, snpQ, snpC, snpD, snpR, phasingR, maxR, max_limit, male) - if state: genotyped_loci += 1 - elif skip_point == 0: print('Locus skipped due to insignificant snps at the level of read split.') - elif skip_point == 1: print('Locus skipped due to less read contribution of Significant snps.') - elif skip_point == 2: print('Locus skipped due to less read contribution in the phased clusters.') - elif skip_point == 6: print(f'Locus {locus_key} skipped due to wide distribution of alleles with one read supporting to it.') # add vcf_fail_writer with new FILTER tag - else: print('Locus skipped due to less number of significant snps based on user\'s parameter.') - else: - if skip_point == 0: - # print('Locus skipped due to minimum number of supporting reads') - vcf_fail_writer(Chrom, locus_key, global_loci_info, ref, out, len(prev_reads), skip_point) - # else: - # print('Locus skipped due to maximum number of supporting reads') - - del global_loci_variations[locus_key] - - return genotyped_loci - -# def flank_adjustment(loci_coords, exact_loci_coords, right_flank_list, left_flank_list): # incorrect code -# # Adjusting flanks based on the nearby repeat region -# for idx,locus in enumerate(loci_coords): -# current_left = left_flank_list[idx] -# current_right = right_flank_list[idx] -# new_left_flanks = [] -# new_right_flanks = [] -# for exact_locus in exact_loci_coords: -# if exact_locus[0]>locus[1]: break -# elif locus[1]>exact_locus[0]: # after skipping locus (or) after came across the current locus, adjust the right flank dist -# new_right_flanks.append(locus[1]-exact_locus[0] if (locus[1]-exact_locus[0]) < current_right else current_right) -# elif locus[0]locus[1]: - break - if e_idx == tr_idx: - continue - if (e_idx>tr_idx) and (locus[1]>exact_locus[0]): # after skipping locus (or) after came across the current locus, adjust the right flank dist - new_right_flanks.append(locus[1]-exact_locus[0] if (locus[1]-exact_locus[0]) < current_right else current_right) - if (e_idx {Chrom} {Start} {End} Total loci = ", total_loci, file=sys.stderr) - - homozygous = False - ambiguous = False - - global_snp_positions = dict() # tracking the encountered SNPs - global_read_variations = {} # tracking the variations on each read - global_loci_variations = {} # tracking the variation for each locus - global_loci_info = {} # saving the information of each loci - - # tracking the loci - global_loci_ends = []; global_loci_keys = [] - global_read_ends = []; global_read_indices = [] - - prev_reads = set() - sorted_global_snp_list = [] - sorted_global_ins_rpos_set = set() - - read_index = 0 - for read in bam.fetch(Chrom, Start[0], End[1]): - - - # skip read with low mapping quality - if read.mapping_quality < mapq_threshold: - continue - - read_chrom = read.reference_name - read_start = read.reference_start - read_end = read.reference_end - - while len(global_loci_ends) > 0 and read_start > global_loci_ends[0]: - - genotyped_loci_count += locus_processor(global_loci_keys, global_loci_ends, global_loci_variations, global_read_variations, global_snp_positions, prev_reads, sorted_global_snp_list, maxR, minR, ref, Chrom, global_loci_info, out, snpQ, snpC, snpD, snpR, phasingR, tbx, flank, sorted_global_ins_rpos_set, log_bool, logger, male) - - - while len(global_read_ends) > 0 and read_start > global_read_ends[0]: - # if the read is beyond the end of the first read that was tracked - if len(global_loci_ends) > 0 and global_read_ends[0] > global_loci_ends[0]: - # if the initial read useful for the first locus being tracked then it is retained - break - else: - - # remove the read information if the current read is beyond the first read and the locus - popped = global_read_ends.pop(0) - rindex = global_read_indices.pop(0) - if rindex in global_read_variations: - for pos in global_read_variations[rindex]['snps']: - if pos in global_snp_positions: - global_snp_positions[pos]['cov'] -= 1 - del_snps = [] - for pos in global_snp_positions: #!! modify this to terminate iteration after some condition - if global_snp_positions[pos]['cov'] == 0: - del_snps.append(pos) - sorted_global_snp_list.remove(pos) - for snp in del_snps: - del global_snp_positions[snp] - del global_read_variations[rindex] - del del_snps - - if rindex in prev_reads: prev_reads.remove(rindex) - del_ins_pos_idx = 0 - list_rpos = sorted(sorted_global_ins_rpos_set) - for i in list_rpos: - del_ins_pos_idx+=1 - if i > popped: break - del list_rpos[:del_ins_pos_idx] - sorted_global_ins_rpos_set = set(list_rpos) - - - - # if the read is beyond the last locus in the bed file the loop stops - if read_start > end_coord: - while len(global_loci_ends) > 0: - genotyped_loci_count += locus_processor(global_loci_keys, global_loci_ends, global_loci_variations, global_read_variations, global_snp_positions, prev_reads, sorted_global_snp_list, maxR, minR, ref, Chrom, global_loci_info, out, snpQ, snpC, snpD, snpR, phasingR, tbx, flank, sorted_global_ins_rpos_set, log_bool, logger, male) - # process the loci left in global_loci_variation - break - - # information locally saved for each read and all the loci it covers - read_loci_variations = {} - # set of homopolymer positions within the reference part that is covered by the read - homopoly_positions = {} - - # repeat loci covered by the read - loci_coords = []; loci_keys = [] - left_flank_list = []; right_flank_list = [] - - for row in tbx.fetch(read_chrom, read_start, read_end): - - # adjust read start and end based on soft and hard clippings - # soft and hard clippings do not consume the reference bases - - row = row.split('\t') - locus_start = int(row[1]); locus_end = int(row[2]); locus_len = locus_end-locus_start - - - if (locus_start>=Start[0]) and (locus_end<=End[1]): - if locus_start==Start[0]: - if locus_end==Start[1]: pass - else: continue - pass - elif locus_start