Mercurial > repos > mheinzl > variant_analyzer2

comparison mut2read.py @ 11:84a1a3f70407 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 Mon, 15 Feb 2021 21:53:24 +0000
parents e18c5293aac7
children d21960b45a6b
comparison
equal deleted inserted replaced
10:e18c5293aac7 11:84a1a3f70407
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 0.2.1 2019-10-27 Gundula Povysil - 14 2.0.0 2020-10-30 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)
66 bam = pysam.AlignmentFile(file2, "rb") 65 bam = pysam.AlignmentFile(file2, "rb")
67 66
68 # get tags 67 # get tags
69 tag_dict = {} 68 tag_dict = {}
70 cvrg_dict = {} 69 cvrg_dict = {}
72 for variant in VCF(file1): 71 for variant in VCF(file1):
73 chrom = variant.CHROM 72 chrom = variant.CHROM
74 stop_pos = variant.start 73 stop_pos = variant.start
75 #chrom_stop_pos = str(chrom) + "#" + str(stop_pos) 74 #chrom_stop_pos = str(chrom) + "#" + str(stop_pos)
76 ref = variant.REF 75 ref = variant.REF
77 if len(variant.ALT) == 0: 76 alt = variant.ALT[0]
78 continue
79 else:
80 alt = variant.ALT[0]
81 print(alt)
82 chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt 77 chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt
83
84
85 dcs_len = [] 78 dcs_len = []
86 if len(ref) == len(alt): 79 if len(ref) == len(alt):
87 for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000): 80 for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000):
88 if pileupcolumn.reference_pos == stop_pos: 81 if pileupcolumn.reference_pos == stop_pos:
89 count_alt = 0 82 count_alt = 0
150 out.write(curr_qual + "\n") 143 out.write(curr_qual + "\n")
151 144
152 145
153 if __name__ == '__main__': 146 if __name__ == '__main__':
154 sys.exit(mut2read(sys.argv)) 147 sys.exit(mut2read(sys.argv))
155