【生信技能樹】fa和fq格式文件的shell小練習(xí)

【生信技能樹】fasta和fastq格式文件的shell小練習(xí)
這個(gè)是有機(jī)結(jié)合生物信息學(xué)的linux和數(shù)據(jù)格式的練習(xí)題：
下載bowtie2軟件后拿到示例數(shù)據(jù)：

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

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

2）輸出所有的reads_1.fq文件中的標(biāo)識(shí)符(即以@開頭的那一行)

$grep -e '^@r' reads_1.fq #解法1
$awk '{if (NR%4==1) print}' reads_1.fq #解法2
T2
  1. 輸出reads_1.fq文件中的 所有序列信息(即每個(gè)序列的第二行)
$awk '{if (NR%4==2) print}' reads_1.fq
T3

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

$awk '{if (NR%4==3) print}' reads_1.fq
T4

5）輸出質(zhì)量值信息(即每個(gè)序列的第四行)

$awk '{if (NR%4==0) print}' reads_1.fq
T5
  1. 計(jì)算reads_1.fq 文件含有N堿基的reads個(gè)數(shù)
$awk '{if (NR%4==2) print}' reads_1.fq | grep 'N' |wc -l
6429
  1. 統(tǒng)計(jì)文件中reads_1.fq文件里面的序列的堿基總數(shù)
$awk '{if (NR%4==2) print}' reads_1.fq |wc -c
1098399

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

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

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

  1. 確定fastQ堿基質(zhì)量值的類型；
  2. 根據(jù)堿基質(zhì)量值類型計(jì)算Q20對(duì)應(yīng)的數(shù)值，轉(zhuǎn)換成ASCII碼；
  3. 統(tǒng)計(jì)該ASCII碼出現(xiàn)次數(shù)。

代碼如下：

$ head reads_1.fq | \
  awk '{if(NR%4==0)printf("%s",$0);}'| \
  od -A n -t u1 | \
  awk 'BEGIN {min=100;max=0;} \
    { for (i=1;i<=NF;i++) \print
      {if ($i>max) max=$i; \
      if ($i<min) min=$i;} \
  } END \
  { if (max<=74 && min<59) print "Phred+33"; \
    else if (max>73 && min>=64) print "Phred+64"; \
    else if (min>=59 && min<64 && max>73) print "Solexa+64"; \
    else print "Unknown score encoding!";
  }'
Phred+33 #確定堿基質(zhì)量值的類型，Q20對(duì)應(yīng)的數(shù)值是20+33=53
$awk 'BEGIN{printf "%c\n",53}'
5 # 53對(duì)應(yīng)的ASCII碼是“5”
$awk '{if (NR%4==0)print}' reads_1.fq | grep -o 5 |wc -l
21369

參考：
【1】一行腳本判斷你的fastq測(cè)序數(shù)據(jù)的質(zhì)量編碼方式
【2】淺談FastQ和FastA格式

10）統(tǒng)計(jì)reads_1.fq 中測(cè)序堿基質(zhì)量值恰好為Q30的個(gè)數(shù)
Q30對(duì)應(yīng)的數(shù)值是30+33=63

$awk 'BEGIN{printf "%c\n",63}'
? # 63對(duì)應(yīng)的ASCII碼是“?”
$awk '{if (NR%4==0)print}' reads_1.fq | grep -o ? |wc -l
21574

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

$awk '{if (NR%4==2) print}' reads_1.fq rea| cut -b 1 | 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)

$awk '{if (NR%4==1||NR%4==2) print}' reads_1.fq | sed 's/@/>/g' > reads_1.fa
T12
  1. 統(tǒng)計(jì)上述reads_1.fa文件中共有多少條序列
$grep -c '>' reads_1.fa
10000

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

$grep -v '>' reads_1.fa | grep -o '[CGcg]' | wc -l
529983

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

$sed 's/N//g' reads_1.fa > reads_1_NoneN.fa
-----驗(yàn)證-----
$cat reads_1.fa | paste - - | grep 'N' | sed 's/\t/\n/g' | head -2
$cat reads_1.fa | paste - - | grep 'N' | sed 's/\t/\n/g' | head -2 | grep 'N'
$head -2 reads_1_NoneN.fa
$head -2 reads_1_NoneN.fa |grep 'N'
r1含有N


在reads_1_NoneN.fa文件中的r1查找N，返回結(jié)果是空

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

$cat reads_1.fa | paste - - | grep -v 'N' | sed 's/\t/\n/g' > reads_1_NN.fa
$head -2 reads_1_NN.fa

在reads_1_NN.fa中r1序列被刪掉了。

  1. 刪除 reads_1.fa文件中的短于65bp的序列
$cat reads_1.fa | paste - - |awk '{if (length($2)>=65) print}' | \
 sed 's/\t/\n/g' > reads_1_morethan65.fa
-----驗(yàn)證-----
$cat reads_1.fa | paste - - |awk '{print(length($2))}' | \ 
 sort | uniq | awk '$2<65{print}' # 查看短于65bp的序列的長(zhǎng)度分布情況
$cat reads_1.fa | paste - - |awk '{print($1,length($2)"\t"$2)}' | sed 's/\t/\n/g' > reads_1_length.fa #在標(biāo)簽后加上序列長(zhǎng)度
$head -2 reads_1_length.fa 
>r1 122
TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG
$cat reads_1_length.fa |paste - - | awk '$2<65{print}' |head -1 
>r9 55  CCCGATGCTTTTTGAAGTTCGCAGAATCGTATGTGTAGANAATTAAACAAANCCT
#查看第一個(gè)短于65bp的序列，是r9
$head -20 reads_1_morethan65.fa | tail -8 #查看reads_1_morethan65.fa文件
>r7
CCCCGCCACCATCCCGCCGGGCNTGTCCATATCGAGCAGAATGCTGTCCACCATCGGATCGCTGGCAGCCTGTTGCAGACGGGCGATAATGCCGTTGTAACCGGTCATCCCCGAGTACGGCTGCAGCGCCCGCGTCCGGCTGA
>r8
NTGAACAGTAAACGTCTGTTGAGCACATCCTTTAATAAGCAGGGCCAGCGCAGTATCNAGTAGCATATTTTTCATGGTGTTATTCCCGATGCTTTTTG
>r10
TTTTCCGGACACAGTTCCGGATGGTCAGCCCGAAGCACATCAGCAACCCGAACAATACCGGCGACAGCCGGAACTGCCGTTCCGGTGTGCAGATTAATGACAGC
>r11
TCGGNCGTCNNTNTGAAGCGGTTATAAATCTGCTCTTTCGCNGTATCCGTANCGATTTCGGTAAGGTAAACCCCG
#查看reads_1_morethan65.fa文件，r9已被刪除

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

$cat reads_1.fa |paste - - | awk '{dels=substr($2,6,length($2)-10);print($1"\n"dels)}'
-----驗(yàn)證-----
$cat reads_1.fa |paste - - | awk '{dels=substr($2,6,length($2)-10);print($1"\n"dels)}' | head -2 | tail -1 > r1_del.fa
$head -2 reads_1.fa |tail -1 > r1.fa
$grep -f r1_del.fa r1.fa
T18.png

19）刪除 reads_1.fa文件中的長(zhǎng)于125bp的序列

$cat reads_1.fa | paste - - |awk '{if (length($2)<=125) print}' | \ 
 sed 's/\t/\n/g' > reads_1_lessthan125.fa
-----驗(yàn)證-----
$cat reads_1_length.fa |paste - - | awk '$2>125{print}' |head -1 
>r2 275 NTTNTGATGCGGGCTTGTGGAGTTCAGCCGATCTGACTTATGTCATTACCTATGAAATGTGAGGACGCTATGCCTGTACCAAATCCTACAATGCCGGTGAAAGGTGCCGGGATCACCCTGTGGGTTTATAAGGGGATCGGTGACCCCTACGCGAATCCGCTTTCAGACGTTGACTGGTCGCGTCTGGCAAAAGTTAAAGACCTGACGCCCGGCGAACTGACCGCTGAGNCCTATGACGACAGCTATCTCGATGATGAAGATGCAGACTGGACTGC
#查看第一個(gè)長(zhǎng)于125bp的序列，是r2
$head reads_1_lessthan125.fa
>r1
TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAAACNGTACGCTGAGGGCAGAAAAAATCGTCGGGGACATTNTAAAGGCGGCGAGCGCGGCTTTTCCG
>r8
NTGAACAGTAAACGTCTGTTGAGCACATCCTTTAATAAGCAGGGCCAGCGCAGTATCNAGTAGCATATTTTTCATGGTGTTATTCCCGATGCTTTTTG
>r9
CCCGATGCTTTTTGAAGTTCGCAGAATCGTATGTGTAGANAATTAAACAAANCCT
>r10
TTTTCCGGACACAGTTCCGGATGGTCAGCCCGAAGCACATCAGCAACCCGAACAATACCGGCGACAGCCGGAACTGCCGTTCCGGTGTGCAGATTAATGACAGC
>r11
TCGGNCGTCNNTNTGAAGCGGTTATAAATCTGCTCTTTCGCNGTATCCGTANCGATTTCGGTAAGGTAAACCCCG
#查看reads_1_lessthan125.fa文件，r2已被刪除

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

$awk '{if(NR%4==0)print}' reads_1.fq | cut -b 1| \
 awk 'BEGIN \
        {for (ii=0; ii<256; ++ii) \
          { ch=sprintf("%c",ii);ascii[ch]=ii;} \
        }\
      {tot+=ascii[$1]}\
      END{print tot/10000-33}'
16.3621
$/mnt/y/biosoft/FastQC/fastqc reads_1.fq

拿reads_1.fq跑了一下fastqc，可以看到第一位堿基的藍(lán)線比16高一點(diǎn)點(diǎn)。答案應(yīng)該沒問題。

FastQC結(jié)果

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容