fasta和fastq格式文件的shell小练习
下载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
fastq文件格式
1)统计reads_1.fq 文件中共有多少条序列信息
fq文件的特点是：每四行代表一条序列。所以我们统计下该文件一共有多少行再除以4就是多少条序列信息啦~
$ wc -l reads_1.fq
paste - - - -的意思是，把四行剪切粘贴成一行
cat reads_1.fq | paste - - - - | wc -l
2）输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)
cat reads_1.fq | paste - - - - | cut -f1
把四行粘成一行，行与行之间默认分隔符为 Tab . 然后 cut -f 取第一个字段便是以@开头的那一行
awk '{if(NR%4==1)print}' reads_1.fq
注意：NR是awk里面的一个内置变量，代表文件行号的意思。
比如我们输入awk '{print NR}' reads_1.fq | head
理解awk内置变量的含义是：指定文件的行号
3) 输出reads_1.fq文件中的所有序列信息(即每个序列的第二行)
cat reads_1.fq | paste - - - - | cut -f2
4）输出以‘+’及其后面的描述信息(即每个序列的第三行)
cat reads_1.fq | paste - - - - | cut -f3
5）输出质量值信息(即每个序列的第四行)
cat reads_1.fq | paste - - - - | cut -f4
6) 计算reads_1.fq 文件含有N碱基的reads个数
我们只需要查看第2列就好
#grep中-c参数：计算符合范本样式的列数。
awk '{if(NR%4==2)print}' reads_1.fq | grep -c N
awk '{if(NR%4==2)print}' reads_1.fq | grep N |wc
# wc（字计数）命令是用来显示文件所包含的行数、字数和字节数。
7) 统计文件中reads_1.fq文件里面的序列的碱基总数
awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN]|wc -l
-o： 是按行输出匹配到的内容
wc -l 是统计行号，这样就是碱基总数啦~
awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + | bc
这种方法的总体思想就是：每条序列都是有长度的，我们把所有序列长度加起来，就是序列的碱基总数。
awk中length参数用法：awk '{if(NR%4==2)print length}' reads_1.fq 这样就会输出每条reads的碱基总数，之后我们把它们全加起来就是碱基总数啦~
理解length参数
要怎么加起来呢？
awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d +paste命令用于合并文件的列。会把每个文件以列对列的方式，一列列地加以合并。
参数
-s，则可以将一个文件中的多行数据合并为一行进行显示。
理解命令：paste -s -d +
加起来之后用命令
bc计算一下。
8）计算reads_1.fq 所有的reads中N碱基的总数
awk '{if(NR%4==2)print}' reads_1.fq | grep -o N |wc
9）统计reads_1.fq 中测序碱基质量值恰好为Q20的个数
质量值在第四列，所以是NR%4==0
Q值在第四列所对应的符号
awk '{if(NR%4==0)print}' reads_1.fq | grep -o 5 |wc
10）统计reads_1.fq 中测序碱基质量值恰好为Q30的个数
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
理解cut的-c参数
所有序列第一位碱基的ATCGN分布情况
12）将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)
看fa/fq的区别
fa有两行/fq有四行
$ cat reads_1.fq | paste - - - - | cut -f1,2|tr '\t' '\n'|tr '@' '>' > reads_1.fa
13) 统计上述reads_1.fa文件中共有多少条序列
$ wc reads_1.fa
14）计算reads_1.fa文件中总的碱基序列的GC数量
fa文件中碱基序列在第二行
$ cat reads_1.fa |grep -o [GC]|wc
grep -o 就是把匹配到的 G 和 C 按行输出。
15）删除 reads_1.fa文件中的每条序列的N碱基
$ cat reads_1.fa |tr -d "N"
tr的-d参数就是删除
理解 tr -d
16）删除 reads_1.fa文件中的含有N碱基的序列
grep -v 取没有匹配的行
$ cat reads_1.fa | paste - - | grep -v N | tr '\t' '\n'
# 这条命令就是，先把两行粘成一行，然后取出所有没有N的序列后再变回去：把两行变成一行。
17) 删除 reads_1.fa文件中的短于65bp的序列
$ cat reads_1.fa | paste - -|awk '{if (length($2)>65) print}'|wc
删除短于65bp的序列也就是说保留大于65bp的序列。
理解 awk 中 length 参数
18） 删除 reads_1.fa文件每条序列的前后五个碱基
cut -c 命令来从文件中抽取任意区间内的字符。例如，cut -c 5-来抽取从位置5开始到行尾的每一个字符;'cut -c -5' 来提取最后5个碱基。
$ cat reads_1.fa | paste - - | cut -f2|cut -c 5- | cut -c -5 ？？？？
# 上面是前5个
19）删除 reads_1.fa文件中的长于125bp的序列
$ cat reads_1.fa | paste - -|awk '{if (length($2)<125) print}'| wc
20）查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值
fq 文件中reads的质量值在第四行
$ awk '{if(NR%4==0)print}' reads_1.fq | cut -c 1 | ????
