这个是有机结合生物信息学的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
做题之前的碎碎念：
本来用conda安装了bowtie2，但是进入文件目录没有发现示例文件夹，找了N久没发现（这是我的第一个疑问o(╥﹏╥)o）。于是按照jimmy老师的方法，用wget下载再解压就有了。之前买了腾讯云99一年的服务器，练手可以，下几个软件可以，但是一牵扯分析数据啊什么的，就内存不不足。于是618活动买了进阶版的云服务器，内存做几个样本的分析是够了。但是昨晚下载bowtie2花了近三个小时，最后让电脑run着,我跑去睡了，今早一看，软件下好了。
拿到example数据，先自己查看一下要联系的数据内容什么的吧。
练习数据：reads_1.fq
1) 统计reads_1.fq 文件中共有多少条序列信息
- 方法一：wc -l
fastq格式总共是4行，wc -l命令是统计有多少行，40000/4=10000条序列
-
方法二：grep
2）输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)
- 方法一：
less reads_1.fa | paste - - - - | cut -f 1
cut -f命令
使用教程实例：
-
方法二：grep （个人本来倾向用grep）
3) 输出reads_1.fq文件中的 所有序列信息(即每个序列的第二行)
- 方法一：
cut -f
- 方法二：
awk
awk命令使用详见：https://man.linuxde.net/awk
4）输出以‘+’及其后面的描述信息(即每个序列的第三行)
和上一题类似
less reads_1.fq | paste - - - - |cut -f 3
or
less -SN reads_1.fq |paste - - - -| awk ' {print $3}'
5）输出质量值信息(即每个序列的第四行)
less reads_1.fq | paste - - - - |cut -f 4
or
less -SN reads_1.fq |paste - - - -| awk ' {print $4}'
6) 计算reads_1.fq 文件含有N碱基的reads个数
7) 统计文件中reads_1.fq文件里面的序列的碱基总数
8）计算reads_1.fq 所有的reads中N碱基的总数
9）统计reads_1.fq 中测序碱基质量值恰好为Q20的个数
如何知道Q20对应的ASCII值，详见
https://www.jianshu.com/p/248308513e2e
根据以上图片知道Q20对应的ASCII为5；Q30为？
因此只要找Q20的个数
10）统计reads_1.fq 中测序碱基质量值恰好为Q30的个数
同Q9
11）统计reads_1.fq 中所有序列的第一位碱基的ATCGNatcg分布情况
12）将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)
fq→fa，需将@替换为>，并只提取前两行
13) 统计上述reads_1.fa文件中共有多少条序列
20000/2=10000
14）计算reads_1.fa文件中总的碱基序列的GC数量
个人更倾向方法三，感觉更简单。grep -o : 只输出文件中匹配到的内容并计数；-e：指定字符串作为查找文件内容的范本样式
15）删除 reads_1.fa文件中的每条序列的N碱基
-
方法一
less reads_1.fa | paste - - | cut -f 2 | sed 's/N//g'
-
方法二：
less -SN reads_1.fa |tr -d 'N'
16）删除 reads_1.fa文件中的含有N碱基的序列
- 方法一
less reads_1.fa | paste - - | cut -f 2 | grep -v 'N'
grep -v 反向输出
17) 删除 reads_1.fa文件中的短于65bp的序列
cat reads_1.fa |paste - - |awk '{if (length($2) >= 65) print $0}'|tr '\t' '\n'
18） 删除 reads_1.fa文件每条序列的前后五个碱基
less reads_1.fa |paste - -|awk ' {print $2}'|sed 's/^.....//'|sed 's/.....$//'
19）删除 reads_1.fa文件中的长于125bp的序列
cat reads_1.fa |paste - - |awk '{if (length($2) <= 125) print $0'
20）查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值
这题 perl -alne 那块不太懂(⊙_⊙)? 看了网上别人的作业答案
less -SN reads_1.fq |paste - - - -|awk ' {print $4}'|cut -c 1|perl -alne '{print ord($_)-33}'|paste -s -d+|bc
ord()函数主要用来返回对应字符的ascii码，chr()主要用来表示ascii码对应的字符他的输入时数字，可以用十进制，也可以用十六进制。
网友评论