fasta和fastq格式文件的shell小练习
2019-08-26 本文已影响0人
生信小学弟
这个是有机结合生物信息学的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
