atac_snvs

atac_snvs/1.0/atac_snvs.py

@@ -0,0 +1,195 @@
+#!/usr/bin/env python
+
+
+
+##### ATTRIBUTION #####
+# MBB 498 Project
+# Author: Rachel LaFrance
+
+
+
+##### ABOUT #####
+# This is a Python script intended to identify SNVs from ATAC-Seq BAM files
+
+
+# Functions:
+
+# extract_barcode : will extract CB_Values (barcodes) from the BAM file
+# process_reads : requires BAM file and reference genome FASTA file, identifies SNVs
+# write_tsv : tsv file(s)
+    # snv_file = all SNVs detected
+    # snv_multi = only multiple (>1) barcoded reads per SNV
+# maf_comparison : filters MAF file by "SNP", then compares the generated tsv file against it to confirm matches
+# find_matching_barcodes : looks for any barcodes that picked up on multiple SNVs for further validation
+
+# Output: Multiple tsv files, but maf_comparison is most informative
+
+
+##### SETUP #####
+import pysam
+from pyfaidx import Fasta 
+from collections import defaultdict
+import pandas as pd
+
+
+
+### input and output files ###
+# Change to appropriate input files before running!
+def main():
+    # the destination of the patient BAM file and reference genome to run process_reads function
+    bam_file = "99-13280_subset.bam" #example directory "data/genome_bams/99-13280_subset.bam"
+    reference_fasta = "reference_genome.fa" # reference genome: https://www.bcgsc.ca/downloads/lcr-modules/genome_fastas/grch38.fa
+    output = process_reads(bam_file, reference_fasta)
+
+    # name the output files
+    snv_file = "mpileup.tsv" # outputs all SNVs detected, used just to ensure the code works and is used in the next function
+    snv_multi = "multiple_barcodes.tsv" # outputs SNVs detected by multiple barcodes - just informational, but could be used in the next function
+    write_tsv(snv_file, snv_multi, output)
+
+    # the destination for patient MAF, desired file chosen from write_tsv function above, and named output
+    maf_comparison("99-13280.maf", "mpileup.tsv", "MAF_comparison.tsv")  
+    # example of patient MAF file, snv_file from function above, and example output file for this function
+
+    # optional! looks for a barcode that picked up on multiple SNVs for further validation
+    maf_comparison_file = "MAF_comparison.tsv" # use the file generated above as the input
+    output_file = "matched_barcodes.tsv" # example name of output file of the matched barcodes
+    find_matching_barcodes(comparison_file, output_file)
+
+
+
+# Extract barcodes from BAM using the CB tag
+def extract_barcodes(read):
+    try:
+        return read.get_tag("CB") 
+    except KeyError:
+        return None #to bypass KeyError
+
+
+
+# Identify any SNVs from barcoded reads
+def process_reads(bam_file, reference_fasta):
+    # open the bam file and reference genome
+    bam = pysam.AlignmentFile(bam_file, "rb")
+    fasta = Fasta(reference_fasta) 
+
+    # use defaultdict - nested dict to store output data:
+    # corresponding to chromosome -> position -> reference_base -> query_base -> set(CB barcodes)
+    output = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(set))))
+
+    for read in bam:
+        # checks if read maps to reference, and doesn't map to secondary sites and split alignments
+        if not read.is_unmapped and not read.is_secondary and not read.is_supplementary:
+            chrom = bam.get_reference_name(read.reference_id) # chromosome name
+            reference_base = fasta[chrom][read.reference_start].seq.upper() # reference base at read start position
+            query_base = read.query_sequence.upper()[0]  # query_base : nucleotide from BAM file
+
+            # identify SNVs : remove same reads, and any "N" bases
+            if query_base != reference_base and query_base != 'N':
+                cb_tag = extract_barcodes(read) # extract the CB barcode 
+                output[chrom][read.reference_start][reference_base][query_base].add(cb_tag) # outputs the chr, location, ref and query base, with CB
+
+    bam.close()
+    fasta.close()
+
+    return output
+
+
+
+# OPTIONS:
+# 1 - snv_file : Output tsv file of all SNVs detected from the BAM
+# 2 - snv_multi : Output tsv file with more than one barcoded read per SNV, narrowed results 
+def write_tsv(snv_file, snv_multi, output):
+
+    with open(snv_file, "w") as output_f:
+        output_f.write("Chromosome\tLocation\tReference_Base\tQuery_Base\tNumber_of_Barcodes\tCB_Values\n") # header
+
+        with open(snv_multi, "w") as s_output:
+            s_output.write("Chromosome\tLocation\tReference_Base\tQuery_Base\tNumber_of_Barcodes\tCB_Values\n")
+
+            # Iterate through the chr -> pos -> ref base _> query in data
+            for chrom, positions in output.items():
+                for position, data in positions.items():
+                    for reference_base, query_bases in data.items():
+                        for query_base, barcodes in query_bases.items():
+                            num_barcodes = len(barcodes) # counts number of barcodes
+                            valid_barcodes = [barcode for barcode in barcodes if barcode is not None]  # to filter out "None" values
+                            if valid_barcodes:  # to only proceed with valid barcodes 
+                                cb_values = ",".join(valid_barcodes) # makes string
+                                output_line = f"{chrom}\t{position}\t{reference_base}\t{query_base}\t{num_barcodes}\t{cb_values}\n"
+
+                                if num_barcodes > 1: # only output more than one barcode
+                                    s_output.write(output_line) # snv_multi
+
+                                output_f.write(output_line) # snv_file
+
+
+
+# Compares all results from snv_file checked against the MAF file, and adds additional corresponding columns from the MAF to the file for more info
+def maf_comparison(maf_file, snv_file, output_file):
+    # using pandas to read MAF file
+    maf_data = pd.read_csv(maf_file, sep='\t', comment='#', header=0, dtype=str)
+
+    # filtering the data based on column 10 and "SNP" values
+    filtered_data = maf_data[maf_data['Variant_Type'] == "SNP"].copy()  # avoid SettingWithCopyWarning
+
+    # convert 'Chromosome' to string in order to avoid DtypeWarning message
+    filtered_data['Chromosome'] = filtered_data['Chromosome'].astype(str)
+
+    # open snv_output file with pandas
+    mpileup_data = pd.read_csv(snv_file, sep='\t', dtype=str)
+    columns_to_include = maf_data.columns.tolist() + mpileup_data.columns[2:].tolist()
+
+    # initialize
+    merged_data = pd.DataFrame(columns=columns_to_include)
+
+    for _, mpileup_row in mpileup_data.iterrows():
+        chromosome = mpileup_row['Chromosome']
+        location = float(mpileup_row['Location'])
+        location_max = location + 1 # to account for the script - location/position off by 1 compared to MAF
+
+        # search MAF data based on conditions - matching with mpileup
+        matching_rows = filtered_data[(filtered_data['Chromosome'] == chromosome) &
+                                      (filtered_data['Start_Position'].astype(float) <= location_max) &
+                                      (filtered_data['End_Position'].astype(float) >= location) &
+                                        (filtered_data['Tumor_Seq_Allele2'] == mpileup_row['Query_Base'])]
+
+        if not matching_rows.empty:
+            # merge the matching rows with additional columns from mpileup data, ignore index to prevent unnecessary tabs and incorrect formatting
+            merged_row = pd.concat([matching_rows.iloc[0], mpileup_row[2:]], axis=0, ignore_index=False)
+            merged_data = pd.concat([merged_data, merged_row.to_frame().T], ignore_index=False)
+
+    # merged data to a new tsv file --> "MAF_comparison"
+    # contains all columns from MAF file, with added barcode information at the end including # of reads and lists the barcodes
+    merged_data.to_csv(output_file, sep='\t', index=False)
+
+
+
+# Iterates through CB_Values to identify if multiple SNVs are picked up from the same barcode
+def find_matching_barcodes(maf_comparison_file, output_file):
+    matching_CB_and_Chr_values = {}
+
+    with open(maf_comparison_file, 'r') as comp_f:
+        next(comp_f)  # to skip header
+        for line in comp_f:
+            fields = line.strip().split('\t')
+            output_CB_value = fields[48]  # get CB_Values (49th column)
+            chromosome = fields[4]  # get Chromsome (5th)
+
+            # Add the row to the list corresponding to its CB_Value and Chromosome
+            key = (output_CB_value, chromosome)
+            if key not in matching_CB_and_Chr_values:
+                matching_CB_and_Chr_values[key] = []
+            matching_CB_and_Chr_values[key].append(fields[:13])  # add columns from the MAF for more information
+
+    with open(output_file, 'w') as matching_f:
+        matching_f.write("CB_Value\tChromosome\tMAF_Info\n")
+        for (CB_value, chromosome), rows in matching_CB_and_Chr_values.items():
+            if len(rows) > 1:
+                for row in rows:
+                    matching_f.write("{}\t{}\t{}\n".format(CB_value, chromosome, '\t'.join(row)))
+
+
+
+if __name__ == "__main__":
+    main()
+

atac_snvs/1.0/env.yaml

@@ -0,0 +1,202 @@
+name: atac_snvs
+
+#copied from kostia
+
+channels:
+  - anaconda
+  - conda-forge
+  - bioconda
+  - python
+  - defaults
+dependencies:
+  - _libgcc_mutex=0.1
+  - _openmp_mutex=4.5
+  - aioeasywebdav=2.4.0
+  - aiohttp=3.8.3
+  - aiosignal=1.2.0
+  - amply=0.1.5
+  - appdirs=1.4.4
+  - async-timeout=4.0.2
+  - attmap=0.13.2
+  - backports=1.1
+  - backports.functools_lru_cache=1.6.4
+  - bcrypt=3.2.2
+  - blas=1.1
+  - boto3=1.24.89
+  - botocore=1.27.89
+  - bottleneck=1.3.5
+  - brotlipy=0.7.0
+  - bzip2=1.0.8
+  - c-ares=1.17.1
+  - ca-certificates=2022.9.24
+  - cachetools=5.2.0
+  - certifi=2022.9.24
+  - cffi=1.14.6
+  - charset-normalizer=2.0.10
+  - coincbc=2.10.5
+  - colorama=0.4.4
+  - commonmark=0.9.1
+  - conda=4.11.0
+  - conda-package-handling=1.7.3
+  - configargparse=1.5.3
+  - connection_pool=0.0.3
+  - cryptography=3.4.7
+  - dataclasses=0.8
+  - datrie=0.8.2
+  - defusedxml=0.7.1
+  - dpath=2.0.6
+  - dropbox=11.34.0
+  - expat=2.2.10
+  - filechunkio=1.8
+  - frozenlist=1.3.1
+  - ftputil=5.0.4
+  - future=0.18.2
+  - git=2.23.0
+  - gitdb=4.0.9
+  - google-api-core=2.10.0
+  - google-api-python-client=2.64.0
+  - google-auth=2.12.0
+  - google-auth-httplib2=0.1.0
+  - google-cloud-core=2.3.2
+  - google-cloud-storage=2.5.0
+  - google-crc32c=1.1.2
+  - google-resumable-media=2.4.0
+  - googleapis-common-protos=1.56.4
+  - grpcio=1.42.0
+  - httplib2=0.20.4
+  - icu=68.1
+  - idna=3.3
+  - importlib-metadata=4.11.4
+  - importlib_resources=5.10.0
+  - iniconfig=1.1.1
+  - jinja2=3.1.2
+  - jmespath=1.0.1
+  - jupyter_core=4.11.1
+  - krb5=1.19.2
+  - ld_impl_linux-64=2.35.1
+  - libarchive=3.5.1
+  - libblas=3.9.0
+  - libcblas=3.9.0
+  - libcrc32c=1.1.1
+  - libcurl=7.78.0
+  - libedit=3.1.20191231
+  - libev=4.33
+  - libffi=3.3
+  - libgcc-ng=12.1.0
+  - libgfortran-ng=12.1.0
+  - libgfortran5=12.1.0
+  - libgomp=12.1.0
+  - libiconv=1.16
+  - liblapack=3.9.0
+  - libnghttp2=1.43.0
+  - libopenblas=0.3.21
+  - libprotobuf=3.19.1
+  - libsodium=1.0.18
+  - libsolv=0.7.19
+  - libssh2=1.9.0
+  - libstdcxx-ng=9.3.0
+  - libxml2=2.9.12
+  - logmuse=0.2.6
+  - lz4-c=1.9.3
+  - lzo=2.10
+  - mamba=0.15.2
+  - markupsafe=2.1.1
+  - multidict=6.0.2
+  - ncurses=6.3
+  - numexpr=2.8.1
+  - numpy-base=1.22.3
+  - oauth2client=4.1.3
+  - openblas=0.3.21
+  - openssl=1.1.1q
+  - packaging=21.3
+  - paramiko=2.11.0
+  - pcre=8.44
+  - peppy=0.35.2
+  - perl=5.26.2
+  - pip=21.2.4
+  - pkgutil-resolve-name=1.3.10
+  - plac=1.3.5
+  - pluggy=1.0.0
+  - ply=3.11
+  - prettytable=3.4.1
+  - protobuf=3.19.1
+  - pulp=2.6.0
+  - py=1.11.0
+  - pyasn1=0.4.8
+  - pyasn1-modules=0.2.8
+  - pycosat=0.6.3
+  - pycparser=2.21
+  - pygments=2.13.0
+  - pynacl=1.5.0
+  - pyopenssl=21.0.0
+  - pysftp=0.2.9
+  - pysocks=1.7.1
+  - pytest=7.1.3
+  - python=3.8.10
+  - python-dateutil=2.8.2
+  - python-fastjsonschema=2.16.2
+  - python-irodsclient=1.1.5
+  - python_abi=3.8
+  - pyu2f=0.1.5
+  - pyyaml=6.0
+  - readline=8.1.2
+  - reproc=14.2.1
+  - reproc-cpp=14.2.1
+  - requests=2.27.1
+  - reretry=0.11.1
+  - rich=12.6.0
+  - rsa=4.9
+  - ruamel_yaml=0.15.80
+  - s3transfer=0.6.0
+  - setuptools=58.0.4
+  - six=1.16.0
+  - slacker=0.14.0
+  - smart_open=6.2.0
+  - snakemake-minimal=7.15.2
+  - sqlite=3.37.0
+  - stone=3.3.1
+  - stopit=1.1.2
+  - tk=8.6.11
+  - tomli=2.0.1
+  - toposort=1.7
+  - tqdm=4.62.3
+  - typing-extensions=4.4.0
+  - typing_extensions=4.4.0
+  - ubiquerg=0.6.2
+  - uritemplate=4.1.1
+  - urllib3=1.26.8
+  - veracitools=0.1.3
+  - wcwidth=0.2.5
+  - wheel=0.37.1
+  - xz=5.2.5
+  - yaml=0.2.5
+  - yarl=1.8.1
+  - yte=1.5.1
+  - zlib=1.2.11
+  - zstd=1.5.0
+  - pip:
+    - attrs==21.4.0
+    - docutils==0.18.1
+    - filelock==3.4.2
+    - gitpython==3.1.26
+    - importlib-resources==5.4.0
+    - ipython-genutils==0.2.0
+    - jsonschema==4.3.3
+    - jupyter-core==4.9.1
+    - nbformat==5.1.3
+    - numpy==1.22.0
+    - pandas==2.1.4
+    - psutil==5.9.0
+    - pyfaidx==0.8.1.1
+    - pyparsing==3.0.6
+    - pyrsistent==0.18.0
+    - pysam==0.22.0
+    - pytz==2021.3
+    - ratelimiter==1.2.0.post0
+    - smart-open==5.2.1
+    - smmap==5.0.0
+    - snakemake==6.13.1
+    - tabulate==0.8.9
+    - traitlets==5.1.1
+    - wrapt==1.13.3
+    - zipp==3.7.0