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

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

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

http://ddrv.cn/a/185067

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
  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
  1. 计算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

  1. 统计文件中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
  1. 统计上述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` 取没有匹配的行

  1. 删除 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字符。

参考资料：
[碱基矿工] 从零开始完整学习全基因组测序数据分析：第3节 数据质控
···bash
echo $PATH |grep -o ":"|wc
cat reads_1.fq |paste - - - - |less -SN
cat reads_1.fq |paste - - - - |cut -f 2|grep -o [ATCGNatcg]|wc

http://ascii.911cha.com/

63 ?

53 5

97 a

107 k
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,332评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,508评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,812评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,607评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,728评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,919评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,071评论 3赞 410
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,802评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,256评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,389评论 4赞 332
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,032评论 3赞 316
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2赞 350

推荐阅读更多精彩内容