☰

查找某个 DNA 碱基序列在文件中出现的次数

tag-icon • tag-icon text-processing scripting bioinformatics
查找某个 DNA 碱基序列在文件中出现的次数

aac作业是编写一个名为“countmatches”的 bash 脚本，该脚本将显示某个序列（例如）在指定文件中出现的次数。该脚本应该至少有两个参数，其中第一个参数必须是包含我们给出的有效 DNA 字符串的文件的路径名。其余参数是仅包含任意顺序的碱基a、c、g和 的字符串。t对于每个有效的参数字符串，它将搜索文件中的 DNA 字符串，并计算该参数字符串在 DNA 字符串（即文件）中出现的非重叠次数。

示例序列和输出是，如果字符串aaccgtttgtaaccggaac位于名为 的文件中dnafile，则脚本应按如下方式工作

$ countmatches dnafile ttt
ttt 1

命令为countmatches dnafile ttt，输出为ttt 1，显示ttt出现一次。

这是我的脚本：

#!/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

这是我得到的错误

[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

答案1

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}

输出：

 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}"

输出：

$ ./countmatches dnafile ttt
ttt, 1

回复@Kusalananda 的评论

上述解决方案很重要不重叠字符串中出现的次数。例如：在字符串“acacaca”中，有两个“aca”不重叠出现，以及三个“aca”重叠出现。为了计数重叠出现次数：

#!/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

匹配ttt序列并报告匹配数很容易：

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

或者，如果序列位于文件中：

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

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

因此，您需要做的就是在 bash 脚本中编写这个想法：

#!/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.

chmod u+x countmatches像这样调用脚本（在使其可执行之后）：

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

答案3

对于文件中行中不重叠的碱基，例如

aaccgtttgtaaccggaac 
acacaca

， 尝试

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

对于重叠的，尝试

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

如果您需要每个文件的计数而不是每行的计数，请将CNTs 相加并在该部分中打印出来END。

相关内容