作业原文:fasta和fastq格式文件的shell小练习 | 生信菜鸟团

原始数据准备

#迅雷下载
#https://github.com/BenLangmead/bowtie2/releases/download/v2.4.1/bowtie2-2.4.1-linux-x86_64.zip
cp ./Desktop/bowtie2-2.4.1-linux-x86_64.zip  ./biosoft/bowtie2
cd ./biosoft/bowtie2
unzip bowtie2-2.4.1-linux-x86_64.zip 
cd ~/biosoft/bowtie2/bowtie2-2.4.1-linux-x86_64/example/reads
ls

Q1

统计reads_1.fq 文件中共有多少条序列信息

grep -c "^@r" reads_1.fq

Q2

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

grep  "^@r" reads_1.fq | more

Q3

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

awk 'NR%4==2' reads_1.fq  #取行数/4的余数为2的行
#或者
sed -n '2~4p' reads_1.fq #第2列开始，每四行取一列
#或者
cat reads_1.fq | paste - - - - | cut -f 2 | less -SN

Q4

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

sed -n '3~4p' reads_1.fq | less -SN
#其它方法同上

Q5

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

sed -n '4~4p' reads_1.fq | less -SN
#其它方法同上

Q6

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

sed -n '2~4p' reads_1.fq | grep -c N

Q7

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

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

Q8

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

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

Q9

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

#20+33 对应的ASCⅡ表的键值为5
sed -n '4~4p' reads_1.fq | grep -o  5 | wc

Q10

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

#30+33 对应的ASCⅡ值表的键值为?
sed -n '4~4p' reads_1.fq | grep -o  ? | wc

Q11

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

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

Q12

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

grep -A 1 "^@r" reads_1.fq| grep -v "^-"  | sed 's/@/>/' > reads_1.fa
#或者
cat reads_1.fq | paste - - - - | cut -f 1,2 | tr "\t" "\n" | sed 's/@/>/' |less -SN

Q13

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

grep -c "^@" reads_1.fa

Q14

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

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

Q15

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

sed "s/N//g" reads_1.fa > tmp.fa

Q16

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

cat reads_1.fa |paste - -  |grep -v N | tr "\t" "\n" | more
#或者
sed 'N;s/\n/_/' reads_1.fa | grep -v N | tr "_" "\n" > tmp2.fa

Q17

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

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

Q18

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

cat reads_1.fa | paste - - | awk '{print $1,substr($2,6,(length($2)-10))}' | tr "\t" "\n" | less -SN

Q19

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

cat reads_1.fa | paste - - | awk 'length($2) <= 125 {print}' | tr "\t" "\n" 

Q20

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

cat reads_1.fq | paste - - - - | cut -f 4 | cut -c 1 | sort | uniq -c > tmp
sed 's/\s//g' tmp | cut -c1-3 >times
echo {33..72} | tr " " "\n" > number
paste -d* times number | paste -s -d + | bc
echo "493621/10000-33" | bc