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

欢迎大家分享自己在完成生信技能树linux作业时的笔记心得。
生信人的 linux 20题
http://www.bio-info-trainee.com/2900.html
fasta和fastq格式文件的shell小练习
http://www.bio-info-trainee.com/3575.html
sam和bam格式文件的shell小练习
http://www.bio-info-trainee.com/3578.html
VCF格式文件的shell小练习
http://www.bio-info-trainee.com/3577.html
根据网上的教程，加上自已的理解，做出小练习。首先打开数据

mkdir   ~/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 文件中共有多少条序列信息

less read.fq |grep "@" |wc -l

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

less read.fq |grep "@" |wc -l

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

sed -n '2~4p'
#awk '{if(NR%4==2)print}'

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

sed -n '3~4p'

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

sed -n '4~4p'

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

sed -n '2~4p' read.fq |grep "N" |wc -l

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

sed -n '2~4p' read.fq |wc |xargs |awk -F" " '{print $3}'

awk '{if(NR%4==2) print length}' read.fq |paste -s -d+ |bc
    paste命令用于合并文件的列。会把每个文件以列对列的方式，一列列地加以合并。
    参数-s，则可以将一个文件中的多行数据合并为一行进行显示。

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

sed -n '2~4p' read.fq |grep -o N |wc -l

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

 sed -n '4~4p' read.fq |grep -o "5"|wc -l

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

 sed -n '4~4p' read.fq |grep -o "?"

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

sed -n '2~4p' read.fq |cut -c1 |sort|uniq -c

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

sed -n '1~4s/^@/>/p;2~4p' read.fq >reads_1.fa
#less -SN reads_1.fq| paste - - - - | cut -f 1,2|tr '\t' '\n'|tr '@' '>'  >> reads_1.fa

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

grep ">" reads_1.fa |wc -l

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

grep -v ">" reads_1.fa |grep -o [GC] |wc -l

在使用代码的过程中，发现GC与[GC]的答案不一样，[GC]是统计“GC”在一起的次数；

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

less reads_1.fa|tr -d "N"|grep N

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

grep -v ">" reads_1.fa |sed 's/N//g' |less

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

less reads_1.fa | paste -d"\t" - -|awk -F "\t" '{if(length($2)>65){print $1"\n"$2}}'

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

less reads_1.fa | paste -d"\t" - -|awk -F "\t" '{print $1"\n"substr($2,6,length($0)-10)}'
##有问题

cat reads_1.fa |paste - -  | tee -a paste.fa|cut -f 2 | sed  s/^.....//g |  sed s/.....$//g |paste paste.fa   -| awk '{print $1 "\t" $3}' |tr '\t' '\n'
    #or
less -SN reads_1.fa |paste - -|awk ' {print $2}'|sed 's/^.....//'|sed 's/.....$//'
##|sed 's/^.....//'  将前五个字符替换为空，sed 's/.....$// 将后五个字符替换为空，^.....代表前五个字符，.....$代表后后五个字符```
    #or
less -SN reads_1.fa |paste - -|awk ' {print $2}'|sed 's/ ^.\{5\}//'|sed 's/ \{5\}.$//'

$cat reads_1.fa |paste - - | awk '{dels=substr($2,6,length($2)-10);print($1"\n"dels)}'
-----验证-----
$cat reads_1.fa |paste - - | awk '{dels=substr($2,6,length($2)-10);print($1"\n"dels)}' | head -2 | tail -1 > r1_del.fa
$head -2 reads_1.fa |tail -1 > r1.fa
$grep -f r1_del.fa r1.fa

19.删除 reads_1.fa文件中的长于125bp的序列

less reads_1.fa | paste -d"\t" - -|awk -F "\t" '{if(length($2)>125){print $1"\n"$2}}'

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

cat reads_1.fq|paste - - - -|cut -f 4|cut -c1|perl -alne '{print ord($_)-33}'| paste -s -d+|bc

less -SN reads_1.fq |paste - - - -|awk ' {print $4}'|cut -c 1|perl -alne '{print ord($_)-33}'|paste -s -d+|BC

less -SN reads_1.fq |paste - - - -|awk ' {print $4}'|cut -c 1|perl -alne '{print ord($_)-33}'|paste -s -d+|bc | awk '{print (($0/10000))}'

$awk '{if(NR%4==0)print}' reads_1.fq | cut -b 1| \
 awk 'BEGIN \
        {for (ii=0; ii<256; ++ii) \
          { ch=sprintf("%c",ii);ascii[ch]=ii;} \
        }\
      {tot+=ascii[$1]}\
      END{print tot/10000-33}'
16.3621