Mercurial > repos > mheinzl > variant_analyzer2

comparison mut2read.py @ 43:d21960b45a6b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Single-Molecule-Genetics/VariantAnalyzerGalaxy/tree/master/tools/variant_analyzer commit ee4a8e6cf290e6c8a4d55f9cd2839d60ab3b11c8
author mheinzl
date Tue, 02 Mar 2021 15:32:41 +0000
parents 84a1a3f70407
children 6ccff403db8a
comparison
equal deleted inserted replaced
42:da224c392a54 43:d21960b45a6b
9 all tags of reads that carry the mutation to a user specified output file. 9 all tags of reads that carry the mutation to a user specified output file.
10 Creates fastq file of reads of tags with mutation. 10 Creates fastq file of reads of tags with mutation.
11 11
12 ======= ========== ================= ================================ 12 ======= ========== ================= ================================
13 Version Date Author Description 13 Version Date Author Description
14 2.0.0 2020-10-30 Gundula Povysil - 14 0.2.1 2019-10-27 Gundula Povysil -
15 ======= ========== ================= ================================ 15 ======= ========== ================= ================================
16 16
17 USAGE: python mut2read.py DCS_Mutations.tabular DCS.bam Aligned_Families.tabular Interesting_Reads.fastq tag_count_dict.json 17 USAGE: python mut2read.py DCS_Mutations.tabular DCS.bam Aligned_Families.tabular Interesting_Reads.fastq tag_count_dict.json
18 """ 18 """
19 19
60 60
61 if os.path.isfile(file3) is False: 61 if os.path.isfile(file3) is False:
62 sys.exit("Error: Could not find '{}'".format(file3)) 62 sys.exit("Error: Could not find '{}'".format(file3))
63 63
64 # read dcs bam file 64 # read dcs bam file
65 # pysam.index(file2)
65 bam = pysam.AlignmentFile(file2, "rb") 66 bam = pysam.AlignmentFile(file2, "rb")
66 67
67 # get tags 68 # get tags
68 tag_dict = {} 69 tag_dict = {}
69 cvrg_dict = {} 70 cvrg_dict = {}
71 for variant in VCF(file1): 72 for variant in VCF(file1):
72 chrom = variant.CHROM 73 chrom = variant.CHROM
73 stop_pos = variant.start 74 stop_pos = variant.start
74 #chrom_stop_pos = str(chrom) + "#" + str(stop_pos) 75 #chrom_stop_pos = str(chrom) + "#" + str(stop_pos)
75 ref = variant.REF 76 ref = variant.REF
76 alt = variant.ALT[0] 77 if len(variant.ALT) == 0:
78 continue
79 else:
80 alt = variant.ALT[0]
81 print(alt)
77 chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt 82 chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt
83
84
78 dcs_len = [] 85 dcs_len = []
79 if len(ref) == len(alt): 86 if len(ref) == len(alt):
80 for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000): 87 for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000):
81 if pileupcolumn.reference_pos == stop_pos: 88 if pileupcolumn.reference_pos == stop_pos:
82 count_alt = 0 89 count_alt = 0
143 out.write(curr_qual + "\n") 150 out.write(curr_qual + "\n")
144 151
145 152
146 if __name__ == '__main__': 153 if __name__ == '__main__':
147 sys.exit(mut2read(sys.argv)) 154 sys.exit(mut2read(sys.argv))
155