fasta和fastq格式文件的shell小練習(xí)

這個(gè)是有機(jī)結(jié)合生物信息學(xué)的linux和數(shù)據(jù)格式的練習(xí)題:
下載bowtie2軟件后拿到示例數(shù)據(jù):

mkdir -p ~/biosoft
cd ~/biosoft
wget -c 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 文件中共有多少條序列信息

grep -c '^@r' reads_1.fq

共10000條序列信息。

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

grep '^@r' reads_1.fq

3)輸出reads_1.fq文件中的 所有序列信息(即每個(gè)序列的第二行)

cat reads_1.fq | paste - - - - | cut -f 2 | less -SN
# awk '{if (NR%4 == 2) print}' reads_1.fq | less -SN 

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

cat reads_1.fq | paste - - - - | cut -f 3 | less -SN

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

cat reads_1.fq | paste - - - - | cut -f 4 | less -SN

6)計(jì)算reads_1.fq 文件含有N堿基的reads個(gè)數(shù)

cat reads_1.fq | paste - - - - | cut -f 2 | grep -n 'N' | wc -l 

有6429條reads中包含了N堿基。

7)統(tǒng)計(jì)文件中reads_1.fq文件里面的序列的堿基總數(shù)

cat reads_1.fq | paste - - - - | cut -f 2 | wc -m

堿基總數(shù)為1098399。

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

cat reads_1.fq | paste - - - - | cut -f 2 | grep -o 'N' | wc -l

N堿基總數(shù)為26001。

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

cat reads_1.fq | paste - - - - | cut -f 4 | grep -o '5' | wc -l
# 20 + 33 ---> 53 ---> 5

Q20的個(gè)數(shù)為21369。

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

cat reads_1.fq | paste - - - - | cut -f 4 | grep -o '?' | wc -l
# 30 + 33 ---> 63 ---> ?

Q30的個(gè)數(shù)為21574。

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

cat reads_1.fq | paste - - - - | cut -f 2 | cut -c 1 | sort | uniq -c

12)將reads_1.fq 轉(zhuǎn)為reads_1.fa文件(即將fastq轉(zhuǎn)化為fasta)

cat reads_1.fq | paste - - - - | cut -f 1,2 | tr '\t' '\n' | tr '@' '>' > reads_1.fa
# awk '{if (NR%4 == 1) {print ">" substr($0, 2)}} {if (NR%4 == 2) print}' reads_1.fq > reads_1.fa

13)統(tǒng)計(jì)上述reads_1.fa文件中共有多少條序列

cat reads_1.fa | grep -v '>' | wc -l
# cat reads_1.fa | paste - - | cut -f 2 | wc -l

共10000條序列信息。

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

grep -v '>' reads_1.fa | grep -o '[CG]' | wc -l 
# cat reads_1.fa | paste - - | awk '{print $2}' | grep -o -e G -e C | wc -l

GC數(shù)量為529983。

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

tr -d 'N' < reads_1.fa
# cat reads_1.fa | tr -d 'N' | grep 'N' -c

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

cat reads_1.fa | paste - - | cut -f 2 | sed '/N/d' | wc -l
# cat reads_1.fa | paste - - | awk '{if ($2 !~ "N") print}' | wc -l

17)刪除 reads_1.fa文件中的短于65bp的序列

cat reads_1.fa | paste - - | cut -f 2 | awk '{if (length($1) >= 65) print}' | wc

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

cat reads_1.fa | paste - - | cut -f 2 | sed 's/^.....//' | sed 's/.....$//' | head

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

cat reads_1.fa | paste - - | cut -f 2 | awk '{if (length($1) <= 125) print}' | wc

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

cat reads_1.fq | paste - - - - | cut -f 4 | cut -c 1 | perl -alne '{print ord($_)-33}' | paste -s -d + | bc 
# 涉及到perl單行命令,還需加強(qiáng)學(xué)習(xí)

參考:
http://www.bio-info-trainee.com/3575.html
http://ascii.911cha.com/

提示

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
最后編輯于
?著作權(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ù)。

相關(guān)閱讀更多精彩內(nèi)容

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