fasta和fastq格式文件的shell小練習

這個是有機結合生物信息學的linux和數(shù)據(jù)格式的練習題： 下載bowtie2軟件后拿到示例數(shù)

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)統(tǒng)計reads_1.fq 文件中共有多少條序列信息

fq文件的特點是：每四行代表一條序列。所以我們統(tǒng)計下該文件一共有多少行再除以4就是多少條序列信息

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

paste - - - -的意思是，把四行剪切粘貼成一行
2）輸出所有的reads_1.fq文件中的標識符(即以@開頭的那一行)

1.awk '{if(NR%4==1)print}' reads_1.fq
2.cat reads_1.fq | paste - - - - | cut -f1
cat reads_1.fq | paste - - - - | cut -f1|head
#@r1
#@r2
#@r3
#@r4
#@r5
#@r6
#@r7
#@r8
#@r9
#@r10
  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個數(shù) 
#grep中-c參數(shù)：計算符合范本樣式的列數(shù)。
$ awk '{if(NR%4==2)print}' reads_1.fq | grep -c N  
6429
$ awk '{if(NR%4==2)print}' reads_1.fq | grep  N |wc  
   6429    6429  782897
#已經讀出的記錄數(shù)，就是行號，從1開始
  1. 統(tǒng)計文件中reads_1.fq文件里面的序列的堿基總數(shù) 
1.awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN]|wc -l
2.awk 'if(NR%4==2)print length' reads_1.fq |paste -s -d + |bc
$ awk '{if(NR%4==2)print length}' reads_1.fq |head
122
275
338
184
138
191
143
98
55
104

8）計算reads_1.fq 所有的reads中N堿基的總數(shù)

awk '{if(NR%4==2)print}' reads_1.fq | grep -o N |wc
 26001   26001   52002

9）統(tǒng)計reads_1.fq 中測序堿基質量值恰好為Q20的個數(shù)

$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o 5 |wc
  21369   21369   42738

10）統(tǒng)計reads_1.fq 中測序堿基質量值恰好為Q30的個數(shù)

$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o ? |wc
  21574   21574   43148

11）統(tǒng)計reads_1.fq 中所有序列的第一位堿基的ATCGNatcg分布情況

$ 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 轉為reads_1.fa文件(即將fastq轉化為fasta)

$  cat reads_1.fq | paste - - - - | cut -f1,2|tr '\t' '\n'|tr '@' '>' > reads_1.fa
  1. 統(tǒng)計上述reads_1.fa文件中共有多少條序列
wc reads_1.fa
20000   20000 1167293 reads_1.fa

14）計算reads_1.fa文件中總的堿基序列的GC數(shù)量

$ cat reads_1.fa |cut -f2|grep -o [GC]|wc
 529983  529983 1059966

15）刪除 reads_1.fa文件中的每條序列的N堿基

$ cat reads_1.fa |tr -d "N"
#tr的-d參數(shù)就是刪除
$ cat reads_1.fa |tr -d "N"|head
>r1
TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTTAAAGGCGGCGAGCGCGGCTTTTCCG
>r2
TTTGATGCGGGCTTGTGGAGTTCAGCCGATCTGACTTATGTCATTACCTATGAAATGTGAGGACGCTATGCCTGTACCAAATCCTACAATGCCGGTGAAAGGTGCCGGGATCACCCTGTGGGTTTATAAGGGGATCGGTGACCCCTACGCGAATCCGCTTTCAGACGTTGACTGGTCGCGTCTGGCAAAAGTTAAAGACCTGACGCCCGGCGAACTGACCGCTGAGCCTATGACGACAGCTATCTCGATGATGAAGATGCAGACTGGACTGC

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

#grep -v 取沒有匹配的行
$ cat reads_1.fa |paste --|grep -v N|tr '\t' '\n'|head
>r1
>r2
>r3
>r4
>r5
GTCAGGAAAGTGGTAAAACTGCAACTCAATTACTGCAATGCCCTCGTAATTAAGTGAATTTACAATATCGTCCTGTTCGGAGGGAAGAACGCGGGATGTTCATTCTTCATCACTTTTAATTGATGTATATGCTCTCTT
>r6
>r7
>r8
>r9
  1. 刪除 reads_1.fa文件中的短于65bp的序列
$ cat reads_1.fa | paste - -|awk '{if (length($2)>65) print}'|wc
   7076   14152  992399
(base) yanke18 15:51:36 ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
$ cat reads_1.fa | paste - -|awk '{if (length($2)>65) print}'|head
>r1 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG
>r2 NTTNTGATGCGGGCTTGTGGAGTTCAGCCGATCTGACTTATGTCATTACCTATGAAATGTGAGGACGCTATGCCTGTACCAAATCCTACAATGCCGGTGAAAGGTGCCGGGATCACCCTGTGGGTTTATAAGGGGATCGGTGACCCCTACGCGAATCCGCTTTCAGACGTTGACTGGTCGCGTCTGGCAAAAGTTAAAGACCTGACGCCCGGCGAACTGACCGCTGAGNCCTATGACGACAGCTATCTCGATGATGAAGATGCAGACTGGACTGC
>r3 ATCGCCCGCAGACACCTTCACGCTGGACTGTTTCGGCTTTTACAGCGTCGCTTCATAATCCTTTTTCGCCGCCGCCATCAGCGTGTTGTAATCCGCCTGCAGGATTTTCCCGTCTTTCNGTGCCTTGNTCAGTTCTTCCTGACGGGCGGTATATTTCTGCAGCGGCGTCTGCAGCCGTTCGTNAGCCTTCTGCGCCTCTTCGGTATATTTCAGCCGTGACGCTTCGGTATCGCTCCGCTGCTGCGCATTTTTCTCCTCTTGAGTCTGCTGCTCAGCCTTCTTTCGGGCGGCTTCAAGCGCAAGACGGGCCTTTTCACGATCATCCCAGTAACGCGCCC
>r4 GGGCCAATGCGCTTACTGATGCGGAATTACGCCGTAAGGCCGCAGATGAGCTTGTCCATATGACTGCGAGAATTAACNGTGGTGAGGCGATCCCTGAACCAGTAAAACAACTTCCTGTCATGGGCGGTAGACCTCTAAATCGTGCACAGGCTCTGGCGAAGATCGCAGAAATCAAAGCTAAGTT
>r5 GTCAGGAAAGTGGTAAAACTGCAACTCAATTACTGCAATGCCCTCGTAATTAAGTGAATTTACAATATCGTCCTGTTCGGAGGGAAGAACGCGGGATGTTCATTCTTCATCACTTTTAATTGATGTATATGCTCTCTT
>r6 AGCGACATTCTTCCTCGGTACATAATCTCCTTTGGCGTTTCCCGATGNCCGTCACGCACATGGNATCCCGTGATGACCTCATTAAAAACACGCTGCAATCCCTCCTCATCTTTGCAGGCGTCCGATTTTTTGCGTTGATTTTTTAATGCAGAATATGCAGTTACCGAGATGTTCCGGTATTTGCAAATCGA
>r7 CCCCGCCACCATCCCGCCGGGCNTGTCCATATCGAGCAGAATGCTGTCCACCATCGGATCGCTGGCAGCCTGTTGCAGACGGGCGATAATGCCGTTGTAACCGGTCATCCCCGAGTACGGCTGCAGCGCCCGCGTCCGGCTGA
>r8 NTGAACAGTAAACGTCTGTTGAGCACATCCTTTAATAAGCAGGGCCAGCGCAGTATCNAGTAGCATATTTTTCATGGTGTTATTCCCGATGCTTTTTG
>r10    TTTTCCGGACACAGTTCCGGATGGTCAGCCCGAAGCACATCAGCAACCCGAACAATACCGGCGACAGCCGGAACTGCCGTTCCGGTGTGCAGATTAATGACAGC
>r11    TCGGNCGTCNNTNTGAAGCGGTTATAAATCTGCTCTTTCGCNGTATCCGTANCGATTTCGGTAAGGTAAACCCCG

18） 刪除 reads_1.fa文件每條序列的前后五個堿基

$ cat reads_1.fa |awk '{if($0~/^>/)print $0;else print substr($0,6,length($0)-10)}'|less -S

19）刪除 reads_1.fa文件中的長于125bp的序列

$ cat reads_1.fa | paste - -|awk '{if (length($2)<125) print}'| wc
   6816   13632  559339
(base) yanke18 16:26:42 ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
$ cat reads_1.fa | paste - -|awk '{if (length($2)<125) print}'| head
>r1 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG
>r8 NTGAACAGTAAACGTCTGTTGAGCACATCCTTTAATAAGCAGGGCCAGCGCAGTATCNAGTAGCATATTTTTCATGGTGTTATTCCCGATGCTTTTTG
>r9 CCCGATGCTTTTTGAAGTTCGCAGAATCGTATGTGTAGANAATTAAACAAANCCT
>r10    TTTTCCGGACACAGTTCCGGATGGTCAGCCCGAAGCACATCAGCAACCCGAACAATACCGGCGACAGCCGGAACTGCCGTTCCGGTGTGCAGATTAATGACAGC
>r11    TCGGNCGTCNNTNTGAAGCGGTTATAAATCTGCTCTTTCGCNGTATCCGTANCGATTTCGGTAAGGTAAACCCCG
>r13    GAAAAACATCGCCGCACAGATGGTTAACTTTGACCGTGAGCAGATGCGTCGGNTCGGCAACAACATGC
>r14    TNANTCAGCAACTGCGTGGACTTCAGGTTGTCCTTCNCTTCGCGCGGGGTATGCCCCTTACTGCTTCCTTTACCCATTCCTCACGCTCCATATATGACAAAACCGCCC
>r17    TTCNNNTAAANGCANTCAGCAACGNTTATGTAAAGAAACAGTAAGATAATACTCAACCCGATGTTTGAGTACGGTCATCATCT
>r18    AGCGCAGTGTCACTGCGCGCCTGTGCACTCTGTGGTGCTGCGGCCAGAATGCGGCGGGCCGTTTTCACGGTCATACCGGG
>r21    CCTTCTCCCATCGACGGACGTCCCACATTGGTGACTTTCACCGTGCGGGTGATCACTTCCTTCGCCGTCACCGCCTT

20）查看reads_1.fq 中每條序列的第一位堿基的質量值的平均值

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

>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
#方法2
>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
最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

友情鏈接更多精彩內容