fasta和fastq格式文件的shell小练习

作者: vicLeo | 来源:发表于2020-01-11 20:33 被阅读0次

作业中参考了：https://www.jianshu.com/p/bc1fe435879c

https://www.jianshu.com/p/d10a85f21889

http://www.huanyujianshe.com/p/41b8c87ab9c5

这个是有机结合生物信息学的linux和数据格式的练习题：
下载bowtie2软件后拿到示例数据：

mkdir -p ~/biosoft
cd ~/biosoft
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip 
unzip bowtie2-2.3.4.3-linux-x86_64.zip 
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads

1)统计reads_1.fq 文件中共有多少条序列信息

wc -l reads_1.fq
cat reads_1.fq | paste - - - - |wc -l

11.png

fq文件的特点是：每四行代表一条序列，所以应计算文件中一共有多少行，再除以4便是序列信息了。

paste - - - -的作用是，把四行剪切粘贴成一行
paste的用法

paste [-s][-d <间隔字符>][--help][--version][文件...]

参数

-d<间隔字符>或--delimiters=<间隔字符>，用指定的间隔字符取代跳格字符

-s或--serial ，串列进行而非平行处理，则可以将一个文件中的多行数据合并为一行进行显示

2）输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)

cat reads_1.fq | paste - - - - |cut -f 1

输出reads_1.fq文件中的所有序列信息(即每个序列的第二行)

cat reads_1.fq | paste - - - - | cut -f2

4）输出以‘+’及其后面的描述信息(即每个序列的第三行)

cat reads_1.fq | paste - - - - | cut -f3

5）输出质量值信息(即每个序列的第四行)

cat reads_1.fq | paste - - - - | cut -f4

计算reads_1.fq 文件含有N碱基的**reads个数

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | cut -f 2| grep "N" |wc -l
6429
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep -c N
6429
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq |grep -c N
6429
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep  N |wc
   6429    6429  782897

less命令
-N 显示每行的行号

-S 行过长时间将超出部分舍弃

awk命令

NR 表示已经读出的记录数，就是行号，从1开始

%4除以4

统计文件中reads_1.fq文件里面的序列的碱基总数

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq |paste - - - - |cut -f2| wc
  10000   10000 1098399
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq |paste - - - - |
cut -f2 | grep -o "[ATCGN]" |wc -l
1088399
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN]| wc -l
1088399
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + | bc
1088399

思路：每条序列都是有长度的，若把所有序列长度加起来，就是序列的碱基总数，length参数会输出每条reads的碱基总数，再把它们全加起来就是碱基总数。

paste命令在此起相加的作用。可用于合并文件的列，会把每个文件以列对列的方式，一列列地加以合并。

-s，可以将一个文件中的多行数据合并为一行进行显示

bc命令用于相加

8）计算reads_1.fq 所有的reads中N碱基的总数

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o N |wc
  26001   26001   52002
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  awk '{if (NR%4==2) print}' reads_1.fq | grep -o N |wc
  26001   26001   52002

9）统计reads_1.fq 中测序碱基质量值恰好为Q20的个数

22.png

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o 5 |wc
  21369   21369   42738
 10）统计**reads_1.fq** 中测序**碱基质量值**恰好为`Q30`的个数  ```
```bash
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o ? |wc
  21574   21574   43148
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq| paste - - - - |grep -o ? |wc
  21574   21574   43148

grep命令
-o 或 --only-matching** : 只显示匹配PATTERN 部分。

11）统计reads_1.fq 中所有序列的第一位碱基的ATCGNatcg分布情况

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | cut -c1| sort | uniq -c
   2184 A
   2203 C
   2219 G
   1141 N
   2253 T
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq| paste - - - - |cut -f2|cut -c1 |sort|uniq -c
   2184 A
   2203 C
   2219 G
   1141 N
   2253 T

cut命令

可以将一串字符作为列来显示，字符字段的记法：

N-：从第N个字节、字符、字段到结尾；
N-M：从第N个字节、字符、字段到第M个（包括M在内）字节、字符、字段；
-M：从第1个字节、字符、字段到第M个（包括M在内）字节、字符、字段。

上面是记法，结合下面选项将摸个范围的字节、字符指定为字段：

-b 表示字节；
-c 表示字符；
-f 表示定义字段。

12）将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)

fastq与fasta的差别

[图片上传中...(44.png-4bb799-1578745528466-0)]

44.png

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  cat reads_1.fq| paste - - - - |cut -f1,2| tr '\t' 'n' | tr '@' '>' > reads_1.fa

统计上述reads_1.fa文件中共有多少条序列

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa| cut -f2 |wc -l
10000
(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ wc reads_1.fa
  10000   10000 1167293 reads_1.fa
nl reads_1.fa| grep 'r' | wc -l10000
10000

nl命令

nl命令在linux系统中用来计算文件中行号。nl 可以将输出的文件内容自动的加上行号！

14）计算reads_1.fa文件中总的碱基序列的GC数量

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  cat reads_1.fa| grep -o [GC] | wc

15）删除 reads_1.fa文件中的每条序列的N碱基

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  cat reads_1.fa| tr -d "N"| head

16）删除 reads_1.fa文件中的含有N碱基的序列

思路：先把两行粘成一行，然后取出所有没有N的序列后再变回去：把两行变成一行。

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  cat reads_1.fa | paste - - | grep -v N | tr '\t' '\n'

grep -v` 取没有匹配的行

删除 reads_1.fa文件中的短于65bp的序列

思路：删除短于65bp的序列 = 保留大于65bp的序列

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa | paste - - | awk '{if (length($2)>65) print}' |wc
   3889    7778  976030

18）删除 reads_1.fa文件每条序列的前后五个碱基

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$  cat reads_1.fa | paste - - | cut -f2 |cut -c 5- |  cut -c -5

cut -c 5- 从第5个碱基开始cut

cut -c -5 从倒数第5个碱基开始cut

以上两个cut的位置在输入时不要互换，保持谁在序列的前面谁先cut

19）删除 reads_1.fa文件中的长于125bp的序列

(base) lchen@VM-0-8-ubuntu:~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa | paste - - | awk '{if (length($2)<125) print}'|head

55.png

图中圈出来的数字为序列号，表现为不连续，证明过滤了长于125bp的序列

20）查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值

awk '{if(NR%4==0)print substr($0,1,1)}' reads_1.fq | awk '
BEGIN{for(i=0;i<256;++i) ord[sprintf("%c",i)]=i}BEGIN{count=0}{ count+=(ord[$1]-33)}END{print count,count/NR}'
163621 16.3621

或
awk '{if(NR%4==0)print substr($0,1,1)}' reads_1.fq | awk 'BEGIN{for(i=0;i<256;++i) ord[sprintf("%c",i)]=i}{print ord[$1]-33}'|paste -s -d+|bc
163621

[Linux C 字符串函数 sprintf()、snprintf() 详解]

字符/Ascii码对照

在C/C++语言中，char也是一种普通的scalable类型，除了字长之外，它与short，int，long这些类型没有本质区别，只不过被大家习惯用来表示字符和字符串而已。

于是，使用“%d”或者“%x”打印一个字符，便能得出它的10进制或16进制的ASCII码；反过来，使用“%c”打印一个整数，便可以看到它所对应的ASCII字符。

http://ascii.911cha.com/

63 ?

53 5

97 a

107 k

网友评论

本文标题：fasta和fastq格式文件的shell小练习

本文链接：https://www.haomeiwen.com/subject/ryvvactx.html

延伸阅读

深度阅读

您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！