資料推薦
- 生信菜鳥團(tuán)的淺談FastQ和FastA格式,以及測序數(shù)據(jù)質(zhì)量控制之FastQC
- 生信技能書論壇的blat簡介與格式解讀
- 還有視頻講解：英語視頻；中文視頻
- 生信技能書系列視頻：https://www.bilibili.com/video/av28813815/?p=12
fasta和fastq格式文件的shell小練習(xí)
1)統(tǒng)計(jì)reads_1.fq 文件中共有多少條序列信息
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ ls
longreads.fq reads_12.fq reads_1.fq reads_2.fq simulate.pl
# 第一種：
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq
...
@r10000
GGTGATGCGCGGCTCCGTGCCGCCAAAGCCGTCCGGCACTGACTNGTCGCAG
+
E<**G2F;';H$%9>*0,;0%---<*9-4B7(5A!4C.C,<".5**$<6,:"
# 第二種：
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | wc
40000 40000 228569
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | wc
10000 40000 2285692
# 我也試過這種：
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq |grep '^@' | wc
10219 10219 93042
# 嗯，還不知道問題出在哪里
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq |grep '^@'
···
@r9966
@(9=@B*;&G<4/F#51*>@B3&0H03@.90-"BHH.#7'*74/.?(&&145G'89#*?:?(!"+8@G02*6B<,#+CE9+?-&67*=1/&4A$:G<:;4965D;;)/B=*?B;'6F//1A#"%7+.1D@=/?93B:A3>.<D%69:/G'6),E4(F(41;'"3C)'?BEC;8$H7A?!5D%3D;-.B'%9>/88>9DEA"H8C6#4"5*63=
@r9967
···
# 嗯，然后及發(fā)現(xiàn)一些奇妙的東西
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | less -SN
2）輸出所有的reads_1.fq文件中的標(biāo)識(shí)符(即以@開頭的那一行)
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | cut -f1
# 或者使用awk
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==1)print}' reads_1.fq
3) 輸出reads_1.fq文件中的 所有序列信息(即每個(gè)序列的第二行)
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | cut -f2
4）輸出以‘+’及其后面的描述信息(即每個(gè)序列的第三行)
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | cut -f3
5）輸出質(zhì)量值信息(即每個(gè)序列的第四行)
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fq | paste - - - - | cut -f4
# -c 計(jì)算符合范本樣式的列數(shù)。
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep -c N
6429
vip39@VM-0-15-ubuntu:~/test/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
7) 統(tǒng)計(jì)文件中reads_1.fq文件里面的序列的堿基總數(shù)
# -o 只輸出文件中匹配到的部分。
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN]|wc
1088399 1088399 2176798
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + |bc
8）計(jì)算reads_1.fq 所有的reads中N堿基的總數(shù)
vip39@VM-0-15-ubuntu:~/test/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）統(tǒng)計(jì)reads_1.fq 中測序堿基質(zhì)量值恰好為Q20的個(gè)數(shù)
-
目前Illumina機(jī)器得到的基本是illumina 1.8方案。
堿基質(zhì)量與對(duì)應(yīng)的ASCII字符
vip39@VM-0-15-ubuntu:~/test/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）統(tǒng)計(jì)reads_1.fq 中測序堿基質(zhì)量值恰好為Q30的個(gè)數(shù)
vip39@VM-0-15-ubuntu:~/test/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
11）統(tǒng)計(jì)reads_1.fq 中所有序列的第一位堿基的ATCGNatcg分布情況
vip39@VM-0-15-ubuntu:~/test/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
12）將reads_1.fq 轉(zhuǎn)為reads_1.fa文件(即將fastq轉(zhuǎn)化為fasta)
vip39@VM-0-15-ubuntu:~/test/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
13) 統(tǒng)計(jì)上述reads_1.fa文件中共有多少條序列
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ wc reads_1.fa
20000 20000 1167293 reads_1.fa
14）計(jì)算reads_1.fa文件中總的堿基序列的GC數(shù)量
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |grep -o G|wc
264740 264740 529480
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |grep -o C|wc
265243 265243 530486
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |grep -o [GC]|wc
529983 529983 1059966
15）刪除 reads_1.fa文件中的每條序列的N堿基
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |tr -d "N"
16）刪除 reads_1.fa文件中的含有N堿基的序列
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |paste - -|grep -v N | wc
3571 7142 340122
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |paste - -|grep -v N | tr '\t' '\n'
17) 刪除 reads_1.fa文件中的短于65bp的序列
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat reads_1.fa |paste - -|awk '{if (length($2)>65) print}'|wc
7076 14152 992399
18） 刪除 reads_1.fa文件每條序列的前后五個(gè)堿基
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ head reads_1.fa|paste - - | cut -f2|cut -c5-
# 上面是前5個(gè)
生信技能樹公益視頻合輯：學(xué)習(xí)順序是linux，r，軟件安裝，geo，小技巧，ngs組學(xué)！
請(qǐng)猛戳下面鏈接：
B站鏈接：https://m.bilibili.com/space/338686099
YouTube鏈接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程師入門最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
學(xué)徒培養(yǎng)：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw
