-1

The assignment is to write a bash script named “countmatches” that will display the number of times a certain sequence, such as aac, appears in a specified file. The script should expect at least two arguments in which the first argument has to be the pathname of a file containing a valid DNA string which we are given. The remaining argument(s) are strings containing only the bases a, c, g, and t in any order.  For each valid argument string, it will search the DNA string in the file and count how many non-overlapping occurrences of that argument string are in the DNA string (i.e., the file).

An example sequence and output would be if the string aaccgtttgtaaccggaac is in a file named dnafile, then the script should work as follows

$ countmatches dnafile ttt
ttt 1

with the command being countmatches dnafile ttt and the output being ttt 1, showing that ttt appears once.

This is my script:

#!/bin/bash
for /data/biocs/b/student.accounts/cs132/data/dna_textfiles
do
        count=$grep -o '[acgt][acgt][acgt]' /data/biocs/b/student.accounts/cs132/data/dna_textfiles | wc -w
        echo {$/data/biocs/b/student.accounts/cs132/data/dna_textfiles} ${count}
done

and this is the error I get

[Osama.Chaudry07@cslab5 assignment3]$ ./countmatches /data/biocs/b/student.accounts/cs132/data/dna_textfiles aac
./countmatches: line 6: '/data/biocs/b/student.accounts/cs132/data/dna_textfiles': not a valid identifier
  • 3
    We're not a script-writing service, but people where will be happy to help you when you hit specific issues with a script you've written. – roaima Oct 14 '18 at 15:45
  • dna_textfiles is a file with nothing but a sequence of letters a, c , g, and t. This is the file for which we have to write a script that will show you how many times a certain sequence such as aac comes up. – Chaudry Osama Oct 14 '18 at 16:04
  • @Goro yes that is what the goal is, to enter any sequence that is present in the large dna_textfiles and have the output be how many times that base appears. I wrote a script for it but it isn't achieving what I want it to and I don't understand where I went wrong. – Chaudry Osama Oct 14 '18 at 17:02
  • @Goro all of the repeats in a sequence. For example in the sequence aaccgtttgtaaccggaac, if I were to input the base aac, it would show that it comes up 3 times. – Chaudry Osama Oct 14 '18 at 17:12
  • 1
    Why are you doing sequence matching using grep and not blast or exonerate or some similar bioinformatics tool? – Kusalananda Oct 14 '18 at 22:28
4
cat dna_textfile 
aaccgtttgtaaccggaac 

#!/bin/bash    
dna_file=/path/to/dna_textfiles
printf "\e[31mNucleotide sequence?:";
read -en 3 userInput
while [[ -z "${userInput}" ]]
do
read -en 3 userInput
done

count=$(grep -o "${userInput}" "${dna_file}" | wc -l)

echo "${userInput}", ${count}

output:

 ttt, 1

#!/bin/bash
#set first and second arguments (dnafile and base respectively)

dir=$1
base=$2

count=$(grep -o ${base} ${dir} | wc -l)

echo "${base}", "${count}"

output:

$ ./countmatches dnafile ttt
ttt, 1

Answer to a comment made by @Kusalananda

The solutions above counts for non-overlapping occurrences within a string. For example: in the string "acacaca", there are two "aca" non-overlapping occurrences, and three "aca" overlapping occurrences. In order to count overlapping occurrences:

#!/bin/bash
#set first and second arguments (sequence and base respectively)  
sequence=$1
base=$2
diff_sequence_base=$((${#sequence} - ${#base} | bc))

for ((i=0; i <= ${diff_sequence_base}; i++)); do
       [ ${sequence:i:${#base}} = $base ] && ((count++))

done
echo $base, $count


$ ./countmatches acacaca aca
aca, 3


$ ./countmatches aaccgtttttaaccggaac ttt
ttt, 3
  • 2
    Consider that sequences may overlap, as the three instances of aca in acacaca. – Kusalananda Oct 14 '18 at 21:33
  • The user did ask: "count how many non-overlapping occurrences of that argument string are in the DNA string". So, grep is fine as coded. There are only two non-overlaping 'aca' in acacaca @Kusalananda – Isaac Oct 14 '18 at 21:53
  • 2
    @Isaac Well, the user should really remove the "bioinformatics" tag if they are trying to do analysis of DNA sequences using grep. This really should be done with some variation of the Smith-Waterman algorithm. – Kusalananda Oct 14 '18 at 22:26
  • @Kusalananda: In the OP’s defense, they never applied the [bioinformatics] tag; Jeff Schaller added it.  But also, from a reality-check perspective, the fact that U&L has a [bioinformatics] tag (with over 100 questions) suggests that the community believes that processing gene text with Unix text tools is appropriate and on-topic.   (By contrast, we are not nearly so receptive to requests to parse JSON and XML without appropriate, specialized tools.) – G-Man Oct 14 '18 at 23:07
  • 1
    @G-Man Processing of bioinformatics text by Unix tools is perfectly OK. Doing sequence matching of DNA with grep is not "bioinformatics" though. It's simply not the correct tool for the job. – Kusalananda Oct 15 '18 at 6:06
3

It is easy to match a ttt sequence and report the count of matches:

$ echo 'aaccgtttgtaaccggaac' | grep -o 'ttt' | wc -l

Or, if the sequence in in a file:

$ echo 'aaccgtttgtaaccggaac'>dnafile
$ grep -o 'ttt' dnafile | wc -l
1

$ grep -o 'aac' dnafile | wc -l
3

So, all you need to do is to code this idea in a bash script:

#!/bin/bash
dnafile=${1-./dnafile}                   # Name of the file to read (arg 1)
shift                                    # Erase arg 1.

for pat; do                              # Process all the other line arguments.
    printf '%s ' "$pat"                  # Print the patern used.
    grep -o "$pat" "$dnafile" | wc -l    # Find the count of matches.
done                                     # done.

Call the script (after making it excutable chmod u+x countmatches) like this:

$ ./countmatches dnafile ttt aac ccgtttg ag
ttt 1
aac 3
ccgtttg 1
ag 0
0

For non-overlapping bases in lines in a file like

aaccgtttgtaaccggaac 
acacaca

, try

awk '{print gsub (base, "&")}' base="ttt" file
1
0

For overlapping ones, try

awk '{while (0 < T=index ($0, base)) {CNT++; $0 = substr($0, T+1)}; print CNT+0;  T = CNT = 0}' base="aca" file
0
3

In case you need the count per file in lieu of per line, sum up the CNTs and print out in the END section.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.