作業(yè)中參考了：http://www.reibang.com/p/bc1fe435879c
http://www.reibang.com/p/d10a85f21889
http://www.huanyujianshe.com/p/41b8c87ab9c5
這個(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 文件中共有多少條序列信息
wc -l reads_1.fq
cat reads_1.fq | paste - - - - |wc -l
fq文件的特點(diǎn)是：每四行代表一條序列，所以應(yīng)計(jì)算文件中一共有多少行遂黍，再除以4便是序列信息了典尾。
paste - - - -的作用是蹦渣，把四行剪切粘貼成一行
paste的用法
paste [-s][-d <間隔字符>][--help][--version][文件...]
參數(shù)
-d<間隔字符>或--delimiters=<間隔字符>哩掺，用指定的間隔字符取代跳格字符
-s或--serial ，串列進(jìn)行而非平行處理徽级，則可以將一個(gè)文件中的多行數(shù)據(jù)合并為一行進(jìn)行顯示
2）輸出所有的reads_1.fq文件中的標(biāo)識符(即以@開頭的那一行)
cat reads_1.fq | paste - - - - |cut -f 1
- 輸出reads_1.fq文件中的 所有序列信息(即每個(gè)序列的第二行)
cat reads_1.fq | paste - - - - | cut -f2
4）輸出以‘+’及其后面的描述信息(即每個(gè)序列的第三行)
cat reads_1.fq | paste - - - - | cut -f3
5）輸出質(zhì)量值信息(即每個(gè)序列的第四行)
cat reads_1.fq | paste - - - - | cut -f4
- 計(jì)算reads_1.fq 文件含有N堿基的**reads個(gè)數(shù)
(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 行過長時(shí)間將超出部分舍棄
NR 表示 已經(jīng)讀出的記錄數(shù)察蹲，就是行號，從1開始
%4除以4
- 統(tǒng)計(jì)文件中reads_1.fq文件里面的序列的堿基總數(shù)
(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
思路：每條序列都是有長度的绒净，若把所有序列長度加起來见咒，就是序列的堿基總數(shù)，length參數(shù)會輸出每條reads的堿基總數(shù)挂疆，再把它們?nèi)悠饋砭褪菈A基總數(shù)改览。
paste命令在此起相加的作用＄脱裕可用于合并文件的列宝当，會把每個(gè)文件以列對列的方式，一列列地加以合并胆萧。
-s庆揩，可以將一個(gè)文件中的多行數(shù)據(jù)合并為一行進(jìn)行顯示
bc命令用于相加
8）計(jì)算reads_1.fq 所有的reads中N堿基的總數(shù)
(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）統(tǒng)計(jì)reads_1.fq 中測序堿基質(zhì)量值恰好為Q20的個(gè)數(shù)
(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）統(tǒng)計(jì)**reads_1.fq** 中測序**堿基質(zhì)量值**恰好為`Q30`的個(gè)數(shù) ```
```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）統(tǒng)計(jì)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
可以將一串字符作為列來顯示鸳碧，字符字段的記法：
N-：從第N個(gè)字節(jié)盾鳞、字符、字段到結(jié)尾瞻离；
N-M：從第N個(gè)字節(jié)腾仅、字符、字段到第M個(gè)（包括M在內(nèi)）字節(jié)套利、字符推励、字段鹤耍；
-M：從第1個(gè)字節(jié)、字符验辞、字段到第M個(gè)（包括M在內(nèi)）字節(jié)稿黄、字符、字段跌造。
上面是記法杆怕，結(jié)合下面選項(xiàng)將摸個(gè)范圍的字節(jié)、字符指定為字段：
-b 表示字節(jié)壳贪；
-c 表示字符陵珍；
-f 表示定義字段。
12）將reads_1.fq 轉(zhuǎn)為reads_1.fa文件(即將fastq轉(zhuǎn)化為fasta)
fastq與fasta的差別
(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
- 統(tǒng)計(jì)上述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命令在linux系統(tǒng)中用來計(jì)算文件中行號违施。nl 可以將輸出的文件內(nèi)容自動(dòng)的加上行號互纯！
14）計(jì)算reads_1.fa文件中總的堿基序列的GC數(shù)量
(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` 取沒有匹配的行
- 刪除 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文件每條序列的前后五個(gè)堿基
(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個(gè)堿基開始cut
cut -c -5 從倒數(shù)第5個(gè)堿基開始cut
以上兩個(gè)cut的位置在輸入時(shí)不要互換留潦，保持誰在序列的前面誰先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
圖中圈出來的數(shù)字為序列號，表現(xiàn)為不連續(xù)辣往，證明過濾了長于125bp的序列
20）查看reads_1.fq 中每條序列的第一位堿基的質(zhì)量值的平均值
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 字符串函數(shù) sprintf()兔院、snprintf() 詳解]
字符/Ascii碼對照
在C/C++語言中，char也是一種普通 的scalable類型站削，除了字長之外秆乳，它與short，int钻哩，long這些類型沒有本質(zhì)區(qū)別，只不過被大家習(xí)慣用來表示字符和字符串而已肛冶。
于是街氢，使用“%d”或者“%x”打印一個(gè)字符睦袖，便能得 出它的10進(jìn)制或16進(jìn)制的ASCII碼珊肃；反過來馅笙，使用“%c”打印一個(gè)整數(shù)，便可以看到它所對應(yīng)的ASCII字符董习。
參考資料：
[堿基礦工] 從零開始完整學(xué)習(xí)全基因組測序數(shù)據(jù)分析：第3節(jié) 數(shù)據(jù)質(zhì)控
···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
