Hi I am very new to this so appologies if this is a simple question - I have been trying to figure this out for days to no avail and my python skills are not quite there yet!
I am trying to extract all N positions from novel bacterial sequences which have been aligned to a member of the same genus. I would like the start and end positions of all N motifs eg. GTCAGNNNNNTGGT
Is there an existing tool / how could I go about creating this in python?
Many thanks.
2 answers
Use SeqKit. SeqKit supports Windows/Linux/Mac OS X.
$ echo -en '>seq\nGTCAGNNNNNTGGT\n' | seqkit locate --ignore-case --only-positive-strand --pattern "N+"
seqID patternName pattern strand start end matched
seq N+ N+ + 6 10 NNNNN
aligned result
$ echo -en '>seq\nGTCAGNNNNNTGGT\n' | seqkit locate --ignore-case --only-positive-strand --pattern "N+" | column -t
seqID patternName pattern strand start end matched
seq N+ N+ + 6 10 NNNNN
or from file and save to file
$ seqkit locate --ignore-case --only-positive-strand --pattern "N+" seqs.fa > result.xls
picard ScatterIntervalsByNs https://broadinstitute.github.io/picard/command-line-overview.html#ScatterIntervalsByNs
Writes an interval list based on splitting a reference by Ns. This tool identifies positions in a reference where the bases are 'no-calls' and writes out an interval-list using the resulting coordinates. This can be used to create an interval list for whole genome sequence (WGS) for e.g. scatter-gather purposes, as an alternative to using fixed-length intervals. The number of contiguous nocalls that can be tolerated before creating a break is adjustable from the command line.
Log in to answer this question.
try regular expression: https://docs.python.org/2/library/re.html
Do you want to search for a specific pattern or for every location in which an 'N' is present?
every location at which a N is present. I have something like 36 different strains and would like to produce a list of N locations in each FASTA file and then compare these lists to find the unique N locations for each strain.
So the output would be the chromosomal locations, right? Sounds like a job that can be done using Biopython. What have you tried?
If the sequence is not long you can do it without software, open the file fasta by the wordpad
That's not very helpful.
thanks but the sequences are > 4 million bases
Just use SeqKit. shenwei356 has even provided a detailed example below.