TODO: edit and sort the following text after vacation
# Input 9,181,515x2 reads
# -- starting point: initial_mutants_rep1.trimmed1_failed_trim
# -- generate bam-file containing not-mapped reads --
#I have a sam-file containing 8586773 alignments from read1. I have 9181515 read2. Can I produce the sorted bam files of 8586773 paired-ends on WA-314_m.fna.
#1115132/2=557566 + 8586773 = 9144339 != 9183535
I have the paired read1 and read2, update the bwa commands
samtools view initial_mutants_rep1.trimmed1_failed_trim | cut -f1 > read1_ids.txt
seqtk subseq initial_mutants_rep1.reads1 read1_ids.txt > filtered_initial_mutants_rep1.reads1.fastq
seqtk subseq initial_mutants_rep1.reads2 read1_ids.txt > filtered_initial_mutants_rep1.reads2.fastq
bwa index ../WA-314_locusid.fasta
bwa mem -M ../WA-314_locusid.fasta filtered_initial_mutants_rep1.reads1.fastq filtered_initial_mutants_rep1.reads2.fastq > initial_mutants_rep1_no_insertion_sites.sam
# Convert SAM to BAM
samtools view -Sb initial_mutants_rep1_no_insertion_sites.sam > initial_mutants_rep1_no_insertion_sites.bam
# Sort the BAM file
samtools sort initial_mutants_rep1_no_insertion_sites.bam -o initial_mutants_rep1_sorted_no_insertion_sites.bam
# Index the sorted BAM file
samtools index initial_mutants_rep1_sorted_no_insertion_sites.bam
#Visualization using tools like IGV to inspect the alignment and distribution of reads.
#Analysis and Visualization: perform various analyses, such as coverage checks, quality assessments.
samtools view -b -f 2 initial_mutants_rep1_sorted_no_insertion_sites.bam > initial_mutants_rep1_properly_paired.bam
bedtools genomecov -ibam initial_mutants_rep1_properly_paired.bam -bg > initial_mutants_rep1_coverage.bedgraph
samtools depth initial_mutants_rep1_properly_paired.bam > initial_mutants_rep1_depth.txt
python3 ~/Scripts/draw_coverages.py
#----> TODO_1 (TOMORROW): draw coverage plots of the reads! senden the results to them!
# TODO_2 (TOMORROW): generate the saturation plots tomorrow!)
samtools faidx WA-314_m.fasta "gi|420257081|ref|NZ_AKKR01000009.1|contig_1_1":31429-31582
Are the following 2 sequences related? for example, subsequences, or reverse complementary?
>gi|420257081|ref|NZ_AKKR01000009.1|contig_1_1:31429-31582
CGTGATCAGGTCTTTTTAGCGACTAAGTTCGGTATTGTGCGTGACCCCACTGATCCCAGG
GTACGCGGTGTGAGCAGCCGTCCTGATTATATCCGCCAGTCGGTTGATGGCAGCCTTAAG
CGCTTGGGTGTTGATGTCATTGATCTTTATTACC
>transposon
CTGTCTCTTATACACATCTCAACCATCATCGATGAATTGTGTCTCAAAATCTCTGATGTTAC
ATTGCACAAGATAAAAATATATCATCATGAACAATAAAACTGTCTGCTTACATAAACAGTAA
TACAAGGGGTGTTATGAGCCATATTCAACGGGAAACGTCTTGCTCGAGGCCGCGATTAAATT
CCAACATGGATGCTGATTTATATGGGTATAAATGGGCTCGCGATAATGTCGGGCAATCAGGT
GCGACAATCTATCGATTGTATGGGAAGCCCGATGCGCCAGAGTTGTTTCTGAAACATGGCAA
AGGTAGCGTTGCCAATGATGTTACAGATGAGATGGTCAGACTAAACTGGCTGACGGAATTTA
TGCCTCTTCCGACCATCAAGCATTTTATCCGTACTCCTGATGATGCATGGTTACTCACCACT
GCGATCCCCGGAAAAACAGCATTCCAGGTATTAGAAGAATATCCTGATTCAGGTGAAAATAT
TGTTGATGCGCTGGCAGTGTTCCTGCGCCGGTTGCATTCGATTCCTGTTTGTAATTGTCCTT
TTAACAGCGATCGCGTATTTCGTCTCGCTCAGGCGCAATCACGAATGAATAACGGTTTGGTT
GATGCGAGTGATTTTGATGACGAGCGTAATGGCTGGCCTGTTGAACAAGTCTGGAAAGAAAT
GCATAAACTTTTGCCATTCTCACCGGATTCAGTCGTCACTCATGGTGATTTCTCACTTGATA
ACCTTATTTTTGACGAGGGGAAATTAATAGGTTGTATTGATGTTGGACGAGTCGGAATCGCA
GACCGATACCAGGATCTTGCCATCCTATGGAACTGCCTCGGTGAGTTTTCTCCTTCATTACA
GAAACGGCTTTTTCAAAAATATGGTATTGATAATCCTGATATGAATAAATTGCAGTTTCATT
TGATGCTCGATGAGTTTTTCTAATCAGAATTGGTTAATTGGTTGTAACACTGGCAGAGCATT
ACGCTGACTTGACGGGACGGCGGCTTTGTTGAATAAATCGAACTTTTGCTGAGTTGAAGGAT
CAGATCACGCATCTTCCCGACAACGCAGACCGTTCCGTGGCAAAGCAAAAGTTCAAAATCAC
CAACTGGTCCACCTACAACAAAGCTCTCATCAACCGTGGCGGGGATCCTCTAGAGTCGACCT
GCAGGCATGCAAGCTTCAGGGTTGAGATGTGTATAAGAGACAG
To determine if the paired reads provide evidence of an insertion site, we need to analyze how the reads align to both the reference genome and the transposon sequence. Typically, for transposon sequencing, we expect to see reads that span the junction between the transposon and the insertion site in the genome.
Given the provided read information and the transposon sequence, here are some key points to consider:
Read Alignment to the Reference Genome:
The read is aligned to positions on the reference genome NZ_AKKR01000009 from 31,438 to 31,543 on the reverse strand.
The CIGAR string 106M52S indicates that 106 bases of the read are mapped to the reference, while the remaining 52 bases are soft clipped (not mapped to the reference).
Mate Read Information:
The mate read starts at position 31,408 on the forward strand of the same reference sequence.
The insert size is -135, suggesting the reads are overlapping or very close to each other.
Supplementary Alignments:
There are supplementary alignments on the reference genome:
From 31,212 to 31,269 on the forward strand (58bp)
From 31,438 to 31,543 on the reverse strand (107bp)
Soft Clipped Sequence:
The 52 soft clipped bases from the read may contain the transposon sequence or part of it, indicating a junction between the transposon and the genome.
To check for evidence of an insertion site, we should:
Extract the soft clipped sequence.
Compare it with the transposon sequence to see if it matches or partially matches the transposon sequence.
Extracting the Soft Clipped Sequence
Given the read length is 158bp and 106 bases are mapped, the remaining 52 bases are soft clipped.
The soft clipped sequence should be extracted from the end of the read.
Check if the Soft Clipped Sequence Matches the Transposon Sequence
You can use sequence alignment tools or a simple script to compare the soft clipped sequence with the transposon sequence.
Let's extract and compare the soft clipped sequence manually.
Manual Extraction and Comparison
Assume the read sequence is as follows:
less
Read sequence (158bp): [106 bases aligned][52 bases soft clipped]
To extract the 52 soft clipped bases:
If the read is provided or can be retrieved from the FASTQ file, you can directly extract the sequence.
Compare the extracted sequence with the start of the transposon sequence.
Example Extraction and Comparison Script
Here is a Python script to extract the soft clipped sequence and compare it with the transposon sequence:
python
# Example read sequence (for illustration, actual sequence needed from FASTQ)
read_sequence = "ACGT... (158bp read sequence here)"
# Extract the 52 soft clipped bases from the end
soft_clipped_sequence = read_sequence[-52:]
# Provided transposon sequence
transposon_sequence = ("CTGTCTCTTATACACATCTCAACCATCATCGATGAATTGTGTCTCAAAATCTCTGATGTTAC"
"ATTGCACAAGATAAAAATATATCATCATGAACAATAAAACTGTCTGCTTACATAAACAGTAA"
"TACAAGGGGTGTTATGAGCCATATTCAACGGGAAACGTCTTGCTCGAGGCCGCGATTAAATT"
"CCAACATGGATGCTGATTTATATGGGTATAAATGGGCTCGCGATAATGTCGGGCAATCAGGT"
"GCGACAATCTATCGATTGTATGGGAAGCCCGATGCGCCAGAGTTGTTTCTGAAACATGGCAA"
"AGGTAGCGTTGCCAATGATGTTACAGATGAGATGGTCAGACTAAACTGGCTGACGGAATTTA"
"TGCCTCTTCCGACCATCAAGCATTTTATCCGTACTCCTGATGATGCATGGTTACTCACCACT"
"GCGATCCCCGGAAAAACAGCATTCCAGGTATTAGAAGAATATCCTGATTCAGGTGAAAATAT"
"TGTTGATGCGCTGGCAGTGTTCCTGCGCCGGTTGCATTCGATTCCTGTTTGTAATTGTCCTT"
"TTAACAGCGATCGCGTATTTCGTCTCGCTCAGGCGCAATCACGAATGAATAACGGTTTGGTT"
"GATGCGAGTGATTTTGATGACGAGCGTAATGGCTGGCCTGTTGAACAAGTCTGGAAAGAAAT"
"GCATAAACTTTTGCCATTCTCACCGGATTCAGTCGTCACTCATGGTGATTTCTCACTTGATA"
"ACCTTATTTTTGACGAGGGGAAATTAATAGGTTGTATTGATGTTGGACGAGTCGGAATCGCA"
"GACCGATACCAGGATCTTGCCATCCTATGGAACTGCCTCGGTGAGTTTTCTCCTTCATTACA"
"GAAACGGCTTTTTCAAAAATATGGTATTGATAATCCTGATATGAATAAATTGCAGTTTCATT"
"TGATGCTCGATGAGTTTTTCTAATCAGAATTGGTTAATTGGTTGTAACACTGGCAGAGCATT"
"ACGCTGACTTGACGGGACGGCGGCTTTGTTGAATAAATCGAACTTTTGCTGAGTTGAAGGAT"
"CAGATCACGCATCTTCCCGACAACGCAGACCGTTCCGTGGCAAAGCAAAAGTTCAAAATCAC"
"CAACTGGTCCACCTACAACAAAGCTCTCATCAACCGTGGCGGGGATCCTCTAGAGTCGACCT"
"GCAGGCATGCAAGCTTCAGGGTTGAGATGTGTATAAGAGACAG")
# Check if the soft clipped sequence matches the start of the transposon sequence
if soft_clipped_sequence in transposon_sequence:
print("The soft clipped sequence matches the transposon sequence, indicating a potential insertion site.")
else:
print("No match found. The soft clipped sequence does not match the transposon sequence.")
# Print the sequences for reference
print("Soft clipped sequence:", soft_clipped_sequence)
print("Transposon sequence start:", transposon_sequence[:52])
Conclusion
To determine if the paired reads provide evidence of an insertion site:
Extract the soft clipped sequence from the read.
Compare the soft clipped sequence with the transposon sequence to see if there is a match.
If there is a match, it indicates a junction between the transposon and the genome, suggesting an insertion site.
The analysis will confirm if the reads are evidence of an insertion site by showing the overlap or junction between the transposon and the genomic sequence.
#--------------
#VH00358:89:AAFC5MTM5:1:2202:39325:27126_:N:0:TTTCTCTA+CTCGACG
I have a lot of reads over 20000 reads in certain positions with soft clips (see the attached image). could you write a script check how many reads from initial_mutants_rep1_no_insertion_sites.sam can confirm the NZ_AKKR01000009:31537 is a insertion site?
python check_insertion_site.py
>transposon
CTGTCTCTTATACACATCTCAACCATCATCGATGAATTGTGTCTCAAAATCTCTGATGTTACATTGCACAAGATAAAAATATATCATCATGAACAATAAAACTGTCTGCTTACATAAACAGTAATACAAGGGGTGTTATGAGCCATATTCAACGGGAAACGTCTTGCTCGAGGCCGCGATTAAATTCCAACATGGATGCTGATTTATATGGGTATAAATGGGCTCGCGATAATGTCGGGCAATCAGGTGCGACAATCTATCGATTGTATGGGAAGCCCGATGCGCCAGAGTTGTTTCTGAAACATGGCAAAGGTAGCGTTGCCAATGATGTTACAGATGAGATGGTCAGACTAAACTGGCTGACGGAATTTATGCCTCTTCCGACCATCAAGCATTTTATCCGTACTCCTGATGATGCATGGTTACTCACCACTGCGATCCCCGGAAAAACAGCATTCCAGGTATTAGAAGAATATCCTGATTCAGGTGAAAATATTGTTGATGCGCTGGCAGTGTTCCTGCGCCGGTTGCATTCGATTCCTGTTTGTAATTGTCCTTTTAACAGCGATCGCGTATTTCGTCTCGCTCAGGCGCAATCACGAATGAATAACGGTTTGGTTGATGCGAGTGATTTTGATGACGAGCGTAATGGCTGGCCTGTTGAACAAGTCTGGAAAGAAATGCATAAACTTTTGCCATTCTCACCGGATTCAGTCGTCACTCATGGTGATTTCTCACTTGATAACCTTATTTTTGACGAGGGGAAATTAATAGGTTGTATTGATGTTGGACGAGTCGGAATCGCAGACCGATACCAGGATCTTGCCATCCTATGGAACTGCCTCGGTGAGTTTTCTCCTTCATTACAGAAACGGCTTTTTCAAAAATATGGTATTGATAATCCTGATATGAATAAATTGCAGTTTCATTTGATGCTCGATGAGTTTTTCTAATCAGAATTGGTTAATTGGTTGTAACACTGGCAGAGCATTACGCTGACTTGACGGGACGGCGGCTTTGTTGAATAAATCGAACTTTTGCTGAGTTGAAGGATCAGATCACGCATCTTCCCGACAACGCAGACCGTTCCGTGGCAAAGCAAAAGTTCAAAATCACCAACTGGTC
|CACCTACAACAAAGCTCTCATCAACCGTGGCGGGGATCCTCTAGAGTCGACCTGCAGGCATGCAAGCTTCAGGGTTGAGATGTGTATAAGAGACAG
add the step to compare the clipped sequences to the transposon sequences. to locate which subsequences of transposon are inserted?
#In "saturation_graph.png," the x-axis represents the input paired reads number, and the y-axis shows the counts for both the template count and TAs hit. It suggests that deeper sequencing will yield more template counts and TAs hits, particularly for the template counts.
samtools flagstat initial_mutants_rep1_sorted_no_insertion_sites.bam
17836098 + 0 in total (QC-passed reads + QC-failed reads)
17173546 + 0 primary
662552 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
17000489 + 0 mapped (95.32% : N/A)
16337937 + 0 primary mapped (95.13% : N/A)
17173546 + 0 paired in sequencing
8586773 + 0 read1
8586773 + 0 read2
15942110 + 0 properly paired (92.83% : N/A)
16234318 + 0 with itself and mate mapped
103619 + 0 singletons (0.60% : N/A)
167124 + 0 with mate mapped to a different chr
163885 + 0 with mate mapped to a different chr (mapQ>=5)
ME Erkennungs Sequenz 5 ́ AGATGTGTATAAGAGACAG 3 ́
ME Sequenz Komplementär: 3 ́ TCTACACATATTCTCTGTC 5 ́
ME Sequenz Rev Komplementär: 5 ́ CTGTCTCTTATACACATCT 3 ́
Kan-2 FP-1 Forward Primer 5 ́ ACCTACAACAAAGCTCTCATCAACC 3 ́
Kan-2 RP-1 Reverse Primer 5 ́ GCAATGTAACATCAGAGATTTTGAG 3 ́
Ampli-Start Senior et al 5’ ACCTACAACAAAGCTCTCATCAACC
3 ́
Sequenzierungs-Strategie (NEBNext Ultra II DNA Modifikation)
NEB_i5_Tn_01
5'-ACACTCTTTCCCTACACGACGCTCTTCCGATCTACCTACAACAAAGCTCTCATCAACC-3'
NEB_i5_Tn_02
5'-ACACTCTTTCCCTACACGACGCTCTTCCGATCTCACCTACAACAAAGCTCTCATCAAC-3'
NEB_i5_Tn_03
5'-ACACTCTTTCCCTACACGACGCTCTTCCGATCTCCTACAACAAAGCTCTCATCAACCG-3'
>> mix 3 forward primer equimolar to create diversity in start of
read 1 (frameshift at the transposon amplification sequence)
NEB_i7_rc
5'-GACTGGAGTTCAGACGTGTGCTCTTCCGATC-3'
>transposon
CTGTCTCTTATACACATCTCAACCATCATCGATGAATTGTGTCTCAAAATCTCTGATGTTACATTGCACAAGATAAAAATATATCATCATGAACAATAAAACTGTCTGCTTACATAAACAGTAATACAAGGGGTGTTATGAGCCATATTCAACGGGAAACGTCTTGCTCGAGGCCGCGATTAAATTCCAACATGGATGCTGATTTATATGGGTATAAATGGGCTCGCGATAATGTCGGGCAATCAGGTGCGACAATCTATCGATTGTATGGGAAGCCCGATGCGCCAGAGTTGTTTCTGAAACATGGCAAAGGTAGCGTTGCCAATGATGTTACAGATGAGATGGTCAGACTAAACTGGCTGACGGAATTTATGCCTCTTCCGACCATCAAGCATTTTATCCGTACTCCTGATGATGCATGGTTACTCACCACTGCGATCCCCGGAAAAACAGCATTCCAGGTATTAGAAGAATATCCTGATTCAGGTGAAAATATTGTTGATGCGCTGGCAGTGTTCCTGCGCCGGTTGCATTCGATTCCTGTTTGTAATTGTCCTTTTAACAGCGATCGCGTATTTCGTCTCGCTCAGGCGCAATCACGAATGAATAACGGTTTGGTTGATGCGAGTGATTTTGATGACGAGCGTAATGGCTGGCCTGTTGAACAAGTCTGGAAAGAAATGCATAAACTTTTGCCATTCTCACCGGATTCAGTCGTCACTCATGGTGATTTCTCACTTGATAACCTTATTTTTGACGAGGGGAAATTAATAGGTTGTATTGATGTTGGACGAGTCGGAATCGCAGACCGATACCAGGATCTTGCCATCCTATGGAACTGCCTCGGTGAGTTTTCTCCTTCATTACAGAAACGGCTTTTTCAAAAATATGGTATTGATAATCCTGATATGAATAAATTGCAGTTTCATTTGATGCTCGATGAGTTTTTCTAATCAGAATTGGTTAATTGGTTGTAACACTGGCAGAGCATTACGCTGACTTGACGGGACGGCGGCTTTGTTGAATAAATCGAACTTTTGCTGAGTTGAAGGATCAGATCACGCATCTTCCCGACAACGCAGACCGTTCCGTGGCAAAGCAAAAGTTCAAAATCACCAACTGGTC
|CACCTACAACAAAGCTCTCATCAACCGTGGCGGGGATCCTCTAGAGTCGACCTGCAGGCATGCAAGCTTCAGGGTTGAGATGTGTATAAGAGACAG
AGAGCTTTGTTGTAGGTG -->
GTGGATGTTGTTTCGAGA
TCTCGAAACAACATCCAC
CACCTACAACAAAGCTCT *
AGAGCTTTGTTGTAGG
