<jupyter_start><jupyter_text>Interactive analysis of .bam file **Part1: Reading the merged-tumour.bam file**<jupyter_code>import pysam
import os

tumor_filename = os.path.join(os.getcwd(), 'Files/merged-tumor.bam')
merged_tumor_file = pysam.AlignmentFile(tumor_filename, "rb")<jupyter_output><empty_output><jupyter_text>**Part 2: Taking out and inspecting the fields of the first read**<jupyter_code>reads = []
for read in merged_tumor_file:
    reads.append(read)

first_read = reads[0]
print(first_read)
print("Query name (QNAME field): " + first_read.query_name)
print("Flag (FLAG field): " + str(first_read.flag))
print("Reference name (RNAME field): " + str(first_read.reference_name))
print("Position (POS field): " + str(first_read.pos))
print("Mapping quality (MAQ field): " + str(first_read.mapping_quality))
print("Cigar string (CIGAR field): " + str(first_read.cigarstring))
print("Next reference id (MRNM field): " + str(first_read.next_reference_id))
print("Next reference position (MPOS field): " + str(first_read.next_reference_start))
print("Insert size (ISIZE field): " + str(first_read.template_length))
print("Query sequence (SEQ field): " + str(first_read.query_sequence))
print("Query quality (QUAL field): " + str(first_read.query_qualities))
print("Tags (TAG field): " + str(first_read.tags))<jupyter_output>Query name (QNAME field): C0HVYACXX120402:7:1207:5722:57044
Flag (FLAG field): 1187
Reference name (RNAME field): 21
Position (POS field): 9483248
Mapping quality (MAQ field): 27
Cigar string (CIGAR field): 76M
Next reference id (MRNM field): 20
Next reference position (MPOS field): 9483381
Insert size (ISIZE field): 209
Query sequence (SEQ field): TTTTCAAACAGTATCTATGCCTGCCAAATGTGAACATATAAAAAAAAACCAGAATGTGCCATTCTGATTTAAACTG
Query quality (QUAL field): array('B', [28, 28, 27, 29, 31, 30, 31, 31, 29, 31, 35, 30, 29, 31, 34, 30, 29, 23, 41, 32, 20, 30, 29, 34, 34, 29, 30, 31, 30, 30, 30, 33, 33, 26, 39, 12, 25, 19, 32, 30, 35, 28, 35, 33, 23, 33, 35, 36, 30, 38, 33, 41, 34, 35, 31, 33, 23, 30, 30, 36, 27, 32, 29, 34, 35, 41, 33, 31, 33, 29, 32, 32, 31, 31, 31, 34])
Tags (TAG field): [('XA', 'GL000217.1,-110754,76M,1;'), ('BD', 'NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN'), ('MD', '76'), ('RG', '1'), ('BI', 'NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN[...]<jupyter_text>Inspecting the flag field<jupyter_code>print(first_read.flag)<jupyter_output>1187<jupyter_text>The value of the flag is 1187 => 0x04A3, which, by the table presented in the lecture means the following: 0x0001 - the read is paired in sequencing; 0x0002 - the read is mapped in a proper pair; 0x0020 - it's a strand of the mate; 0x0080 - it's the second read in a pair and 0x0400 - the read is either a PCR duplicate or an optical duplicate **Part 3: Check out the flag of some reads - taking reads 10, 15 and 50 as examples**<jupyter_code>print("Read 5 flags: " + str(reads[4].flag))
print("Read 15 flags: " + str(reads[14].flag))
print("Read 50 flags: " + str(reads[49].flag))<jupyter_output>Read 5 flags: 99
Read 15 flags: 147
Read 50 flags: 147<jupyter_text>Quick overview of the sample read flags:99 => 0x0063 - the read is: paired in sequencing, mapped in a proper pair, strand of the mate and the first read in a pair 147 => 0x0093 - the read is: paired in sequencing, mapped in a proper pair, strand of the query and the second read in a pair **Part 4: Calculating the number of unmapped reads**<jupyter_code>print("Number of unmapped reads: " + str(len(list(filter(lambda r: r.is_unmapped, reads)))))<jupyter_output>Number of unmapped reads: 17765<jupyter_text>**Part 5: Calculating the total number of reads**<jupyter_code>print("Total number of reads: " + str(len(reads)))<jupyter_output>Total number of reads: 2921629<jupyter_text>**Part 6: Calculating the number of reads with mapping quality = 0**<jupyter_code>print("Number of reads with mapping quality = 0: " + str(len(list(filter(lambda r: r.mapping_quality == 0, reads)))))<jupyter_output>Number of reads with mapping quality = 0: 126628<jupyter_text>**Part 7: Calculating the average mapping quality**<jupyter_code>print("Average mapping quality including reads with mapping quality == 0: " + str(sum(list(map(lambda r: r.mapping_quality, reads)))/len(reads)))<jupyter_output>Average mapping quality including reads with mapping quality == 0: 55.91379158681681<jupyter_text>**Part 8: Calculating the average mapping quality if reads with mapping quality == 0 are taken out**<jupyter_code>reads_filtered = list(filter(lambda r: r.mapping_quality != 0, reads))
print("Number of reads with mapping quality != 0: " + str(len(reads_filtered)))
print("Average mapping quality without reads with mapping quality == 0: " + str(sum(list(map(lambda r: r.mapping_quality, reads_filtered)))/len(reads_filtered)))<jupyter_output>Number of reads with mapping quality != 0: 2795001
Average mapping quality without reads with mapping quality == 0: 58.446975510921106
{'fraction_non_alphanumeric': 0.1278180354267311, 'fraction_numerical': 0.07669082125603864, 'mean_word_length': 6.416417910447761, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '45996', 'license': 'MIT', 'n_tokens_mistral': 2044, 'n_tokens_neox': 1672, 'n_words': 544, 'path': 'GI2022_DZ2.ipynb', 'repo_name': 'batinicm/My-Data-Science-Center', 'size': '4968'}