I am not sure whether tidyverse supports fuzzy search or grouping and I would like to collapse some sequencing data into groups. This is because the sequencer may generate some artificial Ns at the ends, while the entire middle stretch is still identical. I would like to find an efficient way to isolate these sequences and group them (and count them).
For example,
ATCGACACACACACACACACACAATC
, NTCGACACACACACACACACACAATC
, and NNCGACACACACACACACACACAATN
will be grouped into the same group by sharing the CGACACACACACACACACACAAT
. Importantly, CGACACACACACACACACACAAT
is the entire and representative sequence, and it will be used for the group name.
Some other constraints are
- The Ns should only be found at both ends (at either one end, or both ends). Hence, if the Ns are found in the middle of the reads, they should be categorized as separate groups.
- I am only concerning standard DNA (or amino acid) codes. Hence, DNA will use
ATCG
andN
represents the random base while amino acid will useARNDBCEQZGHILKMFPSTWYV
andX
represents the random residue.
Below is a toy dataset and the expected result:
data<- data.frame(DNA = c("ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "NNCGACACACACACACACACACAATN", "ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "ATCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "NTCGACACACACACACACACACAATC", "AGAGTCTCGATCG", "TCTCTGATGCTAAA", "TAGCTAGACTAGCATCGACTACGACT", "TAGCTAGACNNGCATCGACTACGACT", "TAGCTAGACNAGCATCGACTACGACT", "TAGCTAGACTAGCATCGACTACGACT", "TAGCTAGACTAGCATCGACTACGACT", "NAGCTAGACTAGCATCGACTACGACT", "NNCTAGCATCGACTACGACT", "NNCTAGCATCGACTACNACT"))
grouped_data <- data.frame(DNA_group = c("CGACACACACACACACACACAAT", # as mentioned in the description above
"AGAGTCTCGATCG", # noted that the sequence might have differen length, ideally, the suggested code let me define minimum matched length between 2 sequence, but let say 10 now, while reality will be over 100
"TCTCTGATGCTAAA",
"CTAGCATCGACTACGACT", # same as the first example
"TAGCTAGACNNGCATCGACTACGACT", # Ns were found at the middle
"TAGCTAGACNAGCATCGACTACGACT", # Ns were found at the middle
"CTAGCATCGACTACNACT"), # Ns were found at the middle and no matter Ns were also found at the ends, the displayed results should ideally (but optionally to this question's answer) have the Ns at the end trimmed
count = c(11, 1, 1, 5, 1, 1, 1))
I tried to use grep
or gsub
to make a temporary column and count for that column. However, I am not sure if that can accommodate dynamic grep
(as the group name will be different and I am not sure what would be the longest representative sequence for each group).
I prefer to do the job in R, but I will consider other open-source languages if they do the job efficiently. Thanks!