fasta和fastq格式文件的shell小练习

这个是有机结合生物信息学的linux和数据格式的练习题：

下载bowtie2软件后拿到示例数据：

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.统计reads_1.fq 文件中共有多少条序列信息

###第一种方法：
cat reads_1.fq |paste - - - - | wc
>10000   40000 2285692
###第二种方法：
cat reads_1.fq |sed -n '1~4p'|wc -l
>10000

2.输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)

awk '{if(NR%4==1)print}' reads_1.fq 
###下面这种方法有错误，因为中间存在一些以@开头的行
#grep '^@' reads_1.fq |less 
#这种也可以
nl reads_1.fq|  sed -n 1~4p > 1.txt

3.输出reads_1.fq文件中的所有序列信息(即每个序列的第二行)

awk '{if(NR%4==2)print}' reads_1.fq 
#或
nl reads_1.fq|  sed -n 2~4p > 2.txt

4.输出以‘+’及其后面的描述信息(即每个序列的第三行)

awk '{if(NR%4==3)print}' reads_1.fq 
#或
nl reads_1.fq|  sed -n 3~4p > 3.txt

5.输出质量值信息(即每个序列的第四行)

awk '{if(NR%4==0)print}' reads_1.fq > 4.txt
#或
nl reads_1.fq|  sed -n 4~4p > 4.txt

6.计算reads_1.fq 文件含有N碱基的reads个数

nl reads_1.fq|  sed -n 2~4p | sed -n /N/p |wc -l
#或
nl reads_1.fq| sed -n 2~4p |grep 'N' |sort |uniq -c|wc -l

7.统计文件中reads_1.fq文件里面的序列的碱基总数

nl reads_1.fq|  sed -n 2~4p |wc -l

8.计算reads_1.fq 所有的reads中N碱基的总数

nl reads_1.fq|  sed -n 2~4p | sed -n /N/p |uniq -c|wc

9.统计reads_1.fq 中测序碱基质量值恰好为Q20的个数

目前Illumina机器得到的基本是illumina 1.8方案。

image

cat reads_1.fq | paste - - - - | cut -f4 | grep -o "5" | wc 
#或
awk '{if(NR%4==0)print}' reads_1.fq | grep -o '5' |wc

10.统计reads_1.fq 中测序碱基质量值恰好为Q30的个数

#同上一题
cat reads_1.fq | paste - - - - | cut -f4 | grep -o "?" | wc 
awk '{if(NR%4==0)print}' reads_1.fq | grep -o '?' |wc

11.统计reads_1.fq 中所有序列的第一位碱基的ATCGNatcg分布情况

awk '{if(NR%4==2)print}' reads_1.fq | cut -c1|sort |uniq -c

12.将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)

sed 's/^@/>/g' reads_1.fq|sed -n '1~4p;2~4p'>reads_1.fa

13.统计上述reads_1.fa文件中共有多少条序列

awk '{if(NR%4==1)print}' reads_1.fa|wc

14.计算reads_1.fa文件中总的碱基序列的GC数量

cat reads_1.fa |grep -o [GC]|uniq -c|wc -l

15.删除 reads_1.fa文件中的每条序列的N碱基

sed 's/N//g' reads_1.fa 
#或
cat reads_1.fa |tr -d "N"

16.删除 reads_1.fa文件中的含有N碱基的序列

cat reads_1.fa |paste - -|grep -v N | tr '\t' '\n'

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

cat reads_1.fa | paste - - | awk 'length($2)>65{print $0}' | tr '\t' '\n'

18.删除 reads_1.fa文件每条序列的前后五个碱基

less -S read_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}'

20.查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值

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}'
#或
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

最后编辑于：2019.08.26 16:14:07