fasta和fastq格式文件的shell小练习 2019-05-04

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

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

grep -n "^@" reads_1.fq | wc -l
10219
awk '{if(NR%4==1)print}' reads_1.fq| wc -l
10000
cat reads_1.fq | paste - - - - |wc -l
10000 #发现两次输出不同
grep -n "^@" reads_1.fq |head -40
1:@r1
5:@r2
9:@r3
13:@r4
17:@r5
21:@r6
25:@r7
29:@r8
33:@r9
37:@r10
41:@r11
45:@r12
49:@r13
53:@r14
57:@r15
61:@r16
65:@r17
69:@r18
73:@r19
77:@r20
81:@r21
85:@r22
89:@r23
93:@r24
97:@r25
101:@r26
105:@r27
108:@*??&9+&AC!#*D=?&'-D.-38%)$/C/>&A.;<@:D0DE#&AH->9H6:B%!CF/4"8A/<=)0.@&4HH?%B1/-*$+A53
109:@r28
113:@r29
117:@r30
121:@r31
125:@r32
129:@r33
133:@r34
137:@r35
141:@r36
145:@r37
149:@r38  
#108行不正常!
cat reads_1.fq | paste - - - - |head -40 |less -SN
27 @r27    AAGAAGACTTAAAGATTTATCGCACTTGCTCAAANGCNGGATTGATACATATTT
 28 @r28    ATNACTATCAAGGCCGCCTGAGTGCGGTTTTACCGCATACCAATAACGCTTCAC
#r27与r28之间就没有发现异常
cat reads_1.fq |head -120  #看原始数据
@r26
GNAAAGTTANNTAAAAAACATACAGATAACCATCTGCGGTGA
+
""&%%$#$(#!+!%+"!$**%*"+$&&)()!("%*#'&)**!
@r27
AAGAAGACTTAAAGATTTATCGCACTTGCTCAAANGCNGGATTGATACATATTTTGACCCTGAATTATTTCGCTTGAATTTGAAT
+
@*??&9+&AC!#*D=?&'-D.-38%)$/C/>&A.;<@:D0DE#&AH->9H6:B%!CF/4"8A/<=)0.@&4HH?%B1/-*$+A53
@r28
ATNACTATCAAGGCCGCCTGAGTGCGGTTTTACCGCATACCAATAACGCTTCACTCGANGCGT
+
0D6')>@F#%E(=/9D5)@15%?6G2GD,#9$BHDD%742&>9">(62:C/.2@!@:?G@<4.
#发现r27fq格式的最后一行的质量信息以@开头,感觉不太正常,

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

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

3 输出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

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

cat reads_1.fq | paste - - - - |cut -f2|grep "N"|wc -l
6429
awk '{if(NR%4==2)print}' reads_1.fq | grep -c N 

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

awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN]|wc -l
1088399
cat reads_1.fq | paste - - - - |cut -f2 |grep -o [ATCGN] |wc -l
1088399

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

awk '{if(NR%4==2)print}' reads_1.fq | grep -o N|wc -l
26001
cat reads_1.fq | paste - - - - |cut -f2 |grep -o N |wc -l
26001

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

cat reads_1.fq | paste - - - - |cut -f4 |grep -o "5" |wc -l 
21369
awk '{if(NR%4==0)print}' reads_1.fq | grep -o "5"|wc -l
21369
#质量值在第四列,NR%4==0?不懂
```shell
cat reads_1.fq | paste - - - - |cut -f4 |grep -o "?" |wc -l 

[grep “5”及“?”的原因(https://upload-images.jianshu.io/upload_images/15992481-85ee6aa19039c97a.png?imageMogr2/auto-orient/)

10 统计reads_1.fq 中测序碱基质量值恰好为Q30的个数

cat reads_1.fq | paste - - - - |cut -f4 |grep -o "?" |wc -l 
21574
awk '{if(NR%4==0)print}' reads_1.fq | grep -o "?"|wc -l
21574

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

cat reads_1.fq | paste - - - - |cut -f2 |cut -c1 #cut -c1表示按照第一个切第一个字母

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

less -SN reads_1.fq| paste - - - - | cut -f 1,2|tr '\t' '\n'|tr '@' '>' >> reads_1.fa

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

less -SN reads_1.fa |paste - - |wc -l
10000

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

less -SN reads_1.fa |paste - - |cut -f2 |grep -o [CG]|wc -l
529983
awk '{if(NR%2==0)print}' reads_1.fa | grep -o [GC]|wc -l
529983

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

less -SN reads_1.fa  |tr -d N

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

less -SN reads_1.fa|paste - -|sed -e '/N/d' 

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

less -SN reads_1.fa|paste - -|awk '{if (length($2)>65) print}'

18 删除 reads_1.fa文件每条序列的前后五个碱基
看了网上其他同学的命令,看了下结果好像不对,不会做了,,,,

cat reads_1.fa | paste - - | cut -f2 | head -1
TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG
cat reads_1.fa | paste - - | cut -f2 |cut -c -5 |  cut -c 5- | head -1
T

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

less -SN reads_1.fa|paste - -|awk '{if (length($2)>125) print}'

20 查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值
不会了,还得在学啊

less -SN reads_1.fq|paste - - - -|cut -f4 |cut -c 1|head
    原文作者:盼玉
    原文地址: https://www.jianshu.com/p/61438814e210
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞