搜寻序列

1.拿到的原始序列：

2.去除>chrM:

awk 'BEGIN{FS=" "}{print $2}' seq.txt > seq.01.txt

3.在末尾添加NGG：

发现换行符有问题：

cat -A seq.01.txt

^M是windows的换行符，应该去掉：

去掉：

sed -i 's/\r//g' seq.01.txt

cat -A seq.01.txt

添加NGG:

sed 's/$/&NGG/g' seq.01.txt | head

3.使用 seqkit locate 寻找：

F/R链+没有mismatch+全基因组

fasta=/media/shen/6a524d78-97d1-481c-b068-8116a4d007f8/sun/refdata/gencode_GRCm38/raw_fasta/GRCm38.p6.genome.fa

取第一个试下：

seqkit locate -i -p aagcactgaaaatgcttagaNGG $fasta -j 20

all：

for i in `cat seq.fine.txt`

do

echo $i greping.....

seqkit locate -i -p $i $fasta -j 20 -d >> result.txt

done

卧槽太慢了。

换个方法：

sed -i 's/\r//g' seq.txt

sed 's/$/&NGG/g' seq.txt > seq.fa

sed 's/ /\n/g' seq.fa > seq.fine.fa

小程序ref: https://github.com/ekg/fasta-to-fastq

$ chmod u+x fasta_to_fastq.pl

$ perl fasta_to_fastq.pl reads.fasta > my_converted_fasta.fq

# 转为大写：

$ cat seq.fine.fq | tr 'a-z' 'A-Z' > seq.fine.ok.fq

index=/media/shen/6a524d78-97d1-481c-b068-8116a4d007f8/sun/refdata/gencode_GRCm38/01_bowtie2/genome

bowtie2 -x $index -N 0 --end-to-end -a --reorder -U seq.fine.fq -p 10 -S seq.sam

不行...

---

---

应该用本地blast：

安装：