# This script allows choosing a specific NGS results raw data file (.fastq) and filtration of the results in several steps. # The output files would be saved in the same folder as the input fastq file. #Each sample would have a folder containing sequences from the different filtration steps and a folder named "Final" would contain the filtered data that can be used for further analyses. # Make sure to install all required packages before running this script. import tkinter as tk from tkinter import filedialog import os from Bio import SeqIO # Define barcode sequences and sample names. # For barcodes, it is recommended to include the sequence downstream to the barcode for better specificity. barcodes = { "Sample 1": "ATCACGACCGCCTCCACTAGCATATG", "Sample 2": "ACAGTGGTCGCCTCCACTAGCATATG", "Sample 3": "CAGATCCACGCCTCCACTAGCATATG", "Sample 4": "ACAAACGGCGCCTCCACTAGCATATG", "Sample 5": "ACCCAGCACGCCTCCACTAGCATATG", "Sample 6": "AACCCCTCCGCCTCCACTAGCATATG", "Sample 7": "CCCAACCTCGCCTCCACTAGCATATG", "Sample 8": "CACCACACCGCCTCCACTAGCATATG", "Sample 9": "GAAACCCACGCCTCCACTAGCATATG", "Sample 10": "TGTGACCACGCCTCCACTAGCATATG", "Sample 11": "AGGGTCAACGCCTCCACTAGCATATG", "Sample 12": "AGGAGTGGCGCCTCCACTAGCATATG" } # Create a file dialog box to choose the input file root = tk.Tk() root.withdraw() input_file = filedialog.askopenfilename() # Get the directory and filename from the chosen input file input_dir, input_filename = os.path.split(input_file) # Create the output directory path output_dir = os.path.join(input_dir, "DeMultiplexed_" + os.path.splitext(input_filename)[0]) # Create the output directory if it doesn't exist if not os.path.exists(output_dir): os.mkdir(output_dir) # Initialize a dictionary to store output file handles output_files = {bc: open(os.path.join(output_dir, bc + ".fastq"), "w") for bc in barcodes} # Loop through input file and write to appropriate output file based on barcode for record in SeqIO.parse(input_file, "fastq"): seq = str(record.seq) for bc, bc_seq in barcodes.items(): if seq.startswith(bc_seq): SeqIO.write(record, output_files[bc], "fastq") break # Close all output file handles for f in output_files.values(): f.close() from Bio import SeqIO, Seq from collections import defaultdict import re import os import tkinter as tk from tkinter import filedialog # Prompt the user to choose a folder root = tk.Tk() root.withdraw() folder_path = output_dir # Find all the fastq files in the selected folder fastq_files = [f for f in os.listdir(folder_path) if f.endswith('.fastq')] # Check if there are any fastq files in the folder if not fastq_files: print("No fastq files found in the selected folder.") exit(1) # Process each fastq file in the folder for fastq_file in fastq_files: # Get the base name of the fastq file base_name, ext = os.path.splitext(fastq_file) # Create a new directory with the same name as the fastq file output_dir = os.path.join(folder_path, base_name) if not os.path.exists(output_dir): os.mkdir(output_dir) # Create a new directory for the final files output_dir_final = os.path.join(folder_path, 'Final') if not os.path.exists(output_dir_final): os.mkdir(output_dir_final) # Define input and output file paths - rename the files as desired. input_file_1 = os.path.join(folder_path, fastq_file) output_file_1 = os.path.join(output_dir, base_name + '-Filtered.fastq') input_file_2 = os.path.join(output_dir, base_name + '-Filtered.fastq') output_file_2 = os.path.join(output_dir, base_name + '-21.fastq') input_file_3 = os.path.join(output_dir, base_name + '-21.fastq') output_file_3 = os.path.join(output_dir, base_name + '-collapsed.fasta') input_file_4 = os.path.join(output_dir, base_name + '-collapsed.fasta') output_file_4 = os.path.join(output_dir, base_name + '-noStop.fasta') input_file_5 = os.path.join(output_dir, base_name + '-noStop.fasta') output_file_5 = os.path.join(output_dir, base_name + '-noTEV.fasta') input_file_6 = os.path.join(output_dir, base_name + '-noTEV.fasta') output_file_6 = os.path.join(output_dir, base_name + '-noMotif2.fasta') input_file_7 = os.path.join(output_dir, base_name + '-noMotif2.fasta') output_file_7 = os.path.join(output_dir, base_name + '-noMotif2-multiply.fasta') input_file_8 = os.path.join(output_dir, base_name + '-noMotif2-multiply.fasta') output_file_8 = os.path.join(output_dir_final, base_name + '.fasta') from Bio import SeqIO, Seq from collections import defaultdict import re import os import tkinter as tk from tkinter import filedialog # Set the motif to search for. # motif_1 would be a sequence upstream to the substrate region, to include only sequences containing that sequence of interest. # motif_2 can be defined for removal of specific unwanted motifs or sequences from the data. motif_1 = "GTAGCTAGCCCCACTACCGCC" motif_2 = re.compile(r"GGGGGGG") # Function to reverse complement a sequence - should be used in case the data used is from the reverse sequencing. def reverse_complement(seq): return str(Seq.Seq(seq).reverse_complement()) # Step 1: Filter reads by motif - this step make sure to include in the data only sequences that contains a specific sequence num_reads_processed = 0 num_reads_with_motif_1 = 0 with open(output_file_1, "w") as out_handle: for record in SeqIO.parse(input_file_1, "fastq"): num_reads_processed += 1 if motif_1 in str(record.seq): SeqIO.write(record, out_handle, "fastq") num_reads_with_motif_1 += 1 print(f"Processed {num_reads_processed} reads, found {num_reads_with_motif_1} reads with the motif_1") # Step 2: Extract 21 bp sequences after motif and create new reads - this step extracts 21bp (which corresponds to 7 amino acids substrate used) to get a list of dna sequences of the substrates only. input_reads = SeqIO.parse(output_file_1, "fastq") output_reads = [] for read in input_reads: if motif_1 in read.seq: start = read.seq.find(motif_1) + len(motif_1) end = start + 21 seq = str(read.seq[start:end]) new_read = read[start:end] new_read.id = read.id + "_21bp" new_read.description = "" output_reads.append(new_read) SeqIO.write(output_reads, output_file_2, "fastq") # Step 3: Reverse complement the sequences and collapse identical sequences - in this step the 21bp are processed to result in a fasta file with ranked dna sequences based on their count in the data seq_counts = defaultdict(int) with open(output_file_2, "r") as in_file: for record in SeqIO.parse(in_file, "fastq"): rev_comp_seq = reverse_complement(str(record.seq)) seq_counts[rev_comp_seq] += 1 sequences = [(seq, count) for seq, count in seq_counts.items()] sequences.sort(key=lambda x: x[1], reverse=True) with open(output_file_3, "w") as out_file: for i, (seq, count) in enumerate(sequences, 1): SeqIO.write(SeqIO.SeqRecord(Seq.Seq(seq), id=f"{i}-{count}", description=""), out_file, "fasta") # Step 4: Translate DNA sequences to protein sequences and remove stop codons. with open(output_file_3, "r") as in_file, open(output_file_4, "w") as out_file: for record in SeqIO.parse(in_file, "fasta"): protein_seq = record.seq.translate() if "*" not in protein_seq: SeqIO.write(SeqIO.SeqRecord(protein_seq, id=record.id, description=""), out_file, "fasta") # Step 5: Remove TEV substrates - this step removes any sequences that shares >70% similarity the substrate used as the template for library construction. import re from Bio.Seq import Seq from Bio.Align import PairwiseAligner # Open input and output files with open(output_file_4, "r") as input_file, open(output_file_5, "w") as output_file: # Initialize variables sequence = "" header = "" aligner = PairwiseAligner() # Loop through input file for line in input_file: line = line.strip() # If line is a header line if line.startswith(">"): # Write the previous sequence to output file if it doesn't contain the motif if not any(aligner.score(Seq("ENLYFQM"), Seq(sequence[start:end])) / 7 >= 0.7 for start in range(len(sequence) - 6) for end in range(start + 7, len(sequence) + 1)): output_file.write(header + "\n" + sequence + "\n") # Update header and reset sequence header = line sequence = "" # If line is a sequence line else: sequence += line # Write the last sequence to output file if it doesn't contain the motif if not any(aligner.score(Seq("ENLYFQM"), Seq(sequence[start:end])) / 7 >= 0.7 for start in range(len(sequence) - 6) for end in range(start + 7, len(sequence) + 1)): output_file.write(header + "\n" + sequence + "\n") # Step 6: Remove sequences with unwanted motif # Define the motif_2 to search for motif_2 = re.compile(r"GGGGGGG") #GGGGGGG can be replaced with an unwanted motif such as LxxR in yeast systems # Define the name of the input fasta file input_file = input_file_6 # Define the name of the output fasta file output_file = output_file_6 # Open the output file for writing with open(output_file, "w") as f_out: # Iterate over the sequences in the input fasta file for seq_record in SeqIO.parse(input_file, "fasta"): # Flag to keep track of whether the motif_2_2 is present motif_2_present = False # Iterate over all possible 4 amino acid peptide sequences in the protein sequence for i in range(len(seq_record.seq) - 3): peptide = str(seq_record.seq[i:i + 4]) # Check if the motif_2 is present in the peptide if motif_2.search(peptide): motif_2_present = True break # If the motif_2 is present, skip this sequence and move on to the next one if motif_2_present: continue # If the motif_2 is not present, write the sequence to the output file SeqIO.write(seq_record, f_out, "fasta") # Step 7: Multiply (deCollapse) - this step is done to find similar substrate sequences that were counted seperately as a result of synonymous mutations # Define the name of the input fasta file input_file = input_file_7 # Define the name of the output fasta file output_file = output_file_7 with open(input_file, "r") as f: data = f.readlines() sequences = {} for line in data: if line.startswith(">"): seq_id = line.strip()[1:] rank = int(seq_id.split("-")[1].split("_")[0]) sequences[seq_id] = {"rank": rank, "seq": ""} else: sequences[seq_id]["seq"] += line.strip() output_sequences = {} for seq_id in sequences: rank = sequences[seq_id]["rank"] seq = sequences[seq_id]["seq"] for i in range(1, rank + 1): new_seq_id = f"{seq_id.split('-')[0]}-s_{i}" output_sequences[new_seq_id] = seq with open(output_file, "w") as f: for seq_id in output_sequences: seq = output_sequences[seq_id] f.write(f">{seq_id}\n{seq}\n") # CollapseAA sequences - in this step all similar amino acids substrates are being collapsed and ranked according to their abundance in the data. from collections import Counter import operator input_file = input_file_8 output_file = output_file_8 # Read input file and count sequences seq_counts = Counter() with open(input_file, "r") as infile: seq = "" for line in infile: if line.startswith(">"): if seq: seq_counts[seq] += 1 seq = "" else: seq += line.strip() if seq: seq_counts[seq] += 1 # Sort sequences by count, most abundant first sorted_seqs = sorted(seq_counts.items(), key=operator.itemgetter(1), reverse=True) # Write sorted sequences to output file, excluding sequences that appear only once with open(output_file, "w") as outfile: for i, (seq, count) in enumerate(sorted_seqs): if count > 1: # Skip sequences that appear only once outfile.write(">seq{0}_{1}\n{2}\n".format(i + 1, count, seq)) print("Collapsing and filtering sequences complete.") print("File organization complete.")