fasta和fastq格式文件的shell小練習(xí)
這個是有機(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)計reads_1.fq 文件中共有多少條序列信息
2）輸出所有的reads_1.fq文件中的標(biāo)識符(即以@開頭的那一行)
paste - - - -把每4行變成4列
或者awk '{if(NR%4==1)print}' reads_1.fq,NR%4==0就是第4行
awk里有20多個變量，NR是其中一個變量
- 輸出reads_1.fq文件中的 所有序列信息(即每個序列的第二行)
4）輸出以‘+’及其后面的描述信息(即每個序列的第三行)
5）輸出質(zhì)量值信息(即每個序列的第四行)
- 計算reads_1.fq 文件含有N堿基的reads個數(shù)
- 統(tǒng)計文件中reads_1.fq文件里面的序列的堿基總數(shù)
awk '{if(NR%4==2)print}' reads_1.fq|grep -o [ATCGN]| wc
awk '{if(NR%4==2)print}' reads_1.fq|paste -s -d + | bc
8）計算reads_1.fq 所有的reads中
N堿基的總數(shù)
9）統(tǒng)計reads_1.fq 中測序堿基質(zhì)量值恰好為
Q20的個數(shù)
質(zhì)量值在第四列
10）統(tǒng)計reads_1.fq 中測序堿基質(zhì)量值恰好為
Q30的個數(shù)
同上
11）統(tǒng)計reads_1.fq 中所有序列的第一位
堿基的ATCGNatcg分布情況
12）將reads_1.fq 轉(zhuǎn)為reads_1.fa文件(即將fastq轉(zhuǎn)化為fasta)
fa是2行，fq是4行，即只要目前fq的第1、2行
- 統(tǒng)計上述reads_1.fa文件中共有多少條序列
總有10000條
14）計算reads_1.fa文件中總的堿基序列的GC數(shù)量
15）刪除 reads_1.fa文件中的每條序列的N堿基
刪除N堿基tr -d 'N'
grep后發(fā)現(xiàn)沒有
16）刪除 reads_1.fa文件中的含有N堿基的序列
刪除含N序列grep -v N
先將頭文件>r9998和序列變?yōu)橐恍?，然后在將空格轉(zhuǎn)換為換行符
- 刪除 reads_1.fa文件中的短于65bp的序列
刪除短于65的序列，那就是留下大于65的序列，顯示出來
或者
18） 刪除 reads_1.fa文件每條序列的前后五個堿基
$ awk '{if(NR%2==0){print substr($0,6,length($0)-10)}}' reads_1.fa|head -1
19）刪除 reads_1.fa文件中的長于125bp的序列
cat reads_1.fa |paste - - |awk '{if(length($2)<=125){print $1"\n"$2}}'|tail -2
20）查看reads_1.fq 中每條序列的第一位堿基的質(zhì)量值的平均值
第一位堿基：cut -c1
sed和grep主要是對文本進(jìn)行“行”的操作，awk會把每一列都取一個名字，從第一列開始：分別為$1,$2...$n，
