2019年5月25日，一个普通的周六，正在听的歌——北京东路的日子，开始学习生信，写博客。 说明：阅读生信宝典和查阅文章的总结，原文请关注公众号生信宝典，参考的博文都附有链接，仅供参考。

生信宝典

系列教程

关于编程学习的一些思考

知乎专栏：https://zhuanlan.zhihu.com/Data-Analysis

这篇文章讲述两个问题：系统学习还是遇到问题再找答案？是否要写博客。

第一个问题，两种途径都可以，都可以成为大神，但其一，要付出足够多的时间去写代码，修复bug；其二，要多思考、总结，在系统学习时给自己出题，在遇到问题时弄清楚一行行代码的意义；其三，要有兴趣和信仰，对代码的信仰、遇到bug时的信仰、资源的信仰、博客的信仰，相信自己，也相信资源。第二个问题，参考第一个问题的其三。

这篇文章提到几个可能有用的书和文档：

廖雪峰老师的 python 教程R语言，dplyr data.table ggplot2 包的帮助文档stackoverflow 网站几乎能找到所有编程问题的答案

NGS基础——FASTQ格式解释和质量评估

一些复制粘贴，一些注释总结，仅助于自己记忆（记忆力相当差，把能记的都记了）和理解，当然可能有错误，后续再慢慢改正吧 关掉音乐专心学习

FASTQ文件格式和命名

用gzip压缩，一般我们都用双端测序，返回两个FASTQ文件，左端和右端分别命名为_1或R1, _2或R2。 如：sample_name_1_1.fq.gz，sample_name_1_2.fq.gz（第一个下划线后面的数字为重复，第二个下划线后面的数字指定哪一端） 第一行以@开头，后面是reads 的ID 以及其他信息。第二行为read序列。第三行以+开头，一般后面无内容；若有则为序列名字，与第一行相同。第四行为reads质量值。若该碱基测序出错率p_error 为0.001，则Q为30，换算公式为：Q = -10log(p_error)。而测序数据多采用Phred33编码，所以30+33=63，那么63对应的ASCII码为xx，一般碱基质量从0-40，因此ASCII码从（0+33）到（40+33）。以下表显示更为清楚： 为了检测传输的FASTQ文件是否完整，一般会附带一个MD5校验文件

MD5值

#获取MD5值 ct@ehbio:~/$ md5sum test.fq.gz | tee > MD5.txt # tee: 同时写入文件并在屏幕输出结果 + d41d8cd98f00b204e9800998ecf8427e test.fq.gz #输出文字 #验证MD5值 ct@ehbio:~/$ md5sum -c MD5.txt test.fq.gz: 确定 #输出文字

测序质量评估

测序质量可以从两个水平评估：测序reads数 和 测序碱基数。它们都可以用简单的bash命令进行计算，但不是最高效的方式。

ct@ehbio:~/ $ zcat test.fq.gz #zcat: 不解压缩就可以查看文件内容 @ehbio1 ACCGAGCTTTTTTTTTTTTTT + ????????????????????? @ehbio2 ACCGAGCTTTTTTTTTTTTTT + ????????????????????? @ehbio3 ACCGAGCTTTTTTTTTTTTTT + ????????????????????? #获取行数，然后除以4，得到reads数 ct@ehbio:~/ $ echo $(( `zcat test.fq.gz | wc -l` /4 )) 4 #echo: 这里用于显示变量，还可用于-n，输出后不换行，-e，开启转义； #$(( )): 将其他进制转成十进制数显示出来; #``与$()都用来重组命令行，先执行``或()里的命令，将结果替换出来，再组成新的命令行。 #wc: 用于计算文件的字数、行数、列数等，-l 只显示行数。 #获取碱基数 ct@ehbio:~/ $ zcat test.fq.gz | grep -A 1 '^@' | grep '^[ACGTN]' | wc -c 66 #查看test.fq.gz内容，从中获取第一列开头为@的行及其之后的内容，从中获取开头为ACGTN的行，计算字符数。 #grep -A<显示列数>或--after-context=<显示列数> 除了显示符合范本样式的那一列之外，并显示该列之后的内容。

获得reads数和碱基数后，可以用表格展示。但为了更直观一般用柱状图展示。样品特别多时用直方图展示，优先关注是否存在测序量异常高或低的样品。

测序质量

一般用FastQC软件来评估测序质量，它是一个基于Java的软件包（并不清楚Java是啥），下载后放到环境变量就可使用。

运行方式：fastqc ehbio_1_1.fq.gz ehbio_1_2.fq.gz生成结果和解释如表所示： 左图显示每个碱基的中位质量得分(箱线图中间的蓝线)都比较高，而右图的的碱基质量得分变化较大，测序质量逐渐下降，测序质量差。可能需要进行一定的预处理比如移除低质量碱基。 图横坐标表示平均GC含量，纵坐标表示reads数。左图显示每个read的GC分布(红线)与理论分布(蓝线)相契合，GC含量均一。右图出现了GC含量双峰，表示测序样品可能存在特定的序列污染如混入了引物二聚体，当这一指标异常并且导致后期的序列比对率很低时，需要引起注意。 接头序列含量图。横坐标表示read中碱基的位置，纵坐标表示接头序列含量。从图中可以得知read的前17 bp是干净的，之后的位置出现了比例不等的接头序列，需要去除。另一方面也反映了这个库的质量很差，有可能需要重建。 当样品数目很多时，会得到更多的质量评估图，结果很不直观。因此可以用下面所示的所有样品测序质量的散点图来表示。通过GC含量(横坐标)和碱基测序质量(纵坐标)来判断样品的测序质量。

两个值都是FAIL的样品，需要重点关注。

附grep命令详解：https://www.cnblogs.com/wangcp-2014/p/5146335.html 贴来一些 1.主要参数 [options]主要参数：

-a或–text 不要忽略二进制的数据。 -A<显示列数>或–after-context=<显示列数> 除了显示符合范本样式的那一列之外，并显示该列之后的内容。 -b或–byte-offset 在显示符合范本样式的那一列之前，标示出该列第一个字符的位编号。 -B<显示列数>或–before-context=<显示列数> 除了显示符合范本样式的那一列之外，并显示该列之前的内容。 -c或–count 计算符合范本样式的列数。 -C<显示列数>或–context=<显示列数>或-<显示列数> 除了显示符合范本样式的那一列之外，并显示该列之前后的内容。 -d<进行动作>或–directories=<进行动作> 当指定要查找的是目录而非文件时，必须使用这项参数，否则grep指令将回报信息并停止动作。 -e<范本样式>或–regexp=<范本样式> 指定字符串做为查找文件内容的范本样式。 -E或–extended-regexp 将范本样式为延伸的普通表示法来使用。 -f<范本文件>或–file=<范本文件> 指定范本文件，其内容含有一个或多个范本样式，让grep查找符合范本条件的文件内容，格式为每列一个范本样式。 -F或–fixed-regexp 将范本样式视为固定字符串的列表。 -G或–basic-regexp 将范本样式视为普通的表示法来使用。 -h或–no-filename 在显示符合范本样式的那一列之前，不标示该列所属的文件名称。 -H或–with-filename 在显示符合范本样式的那一列之前，表示该列所属的文件名称。 -i或–ignore-case 忽略字符大小写的差别。 -l或–file-with-matches 列出文件内容符合指定的范本样式的文件名称。 -L或–files-without-match 列出文件内容不符合指定的范本样式的文件名称。 -n或–line-number 在显示符合范本样式的那一列之前，标示出该列的列数编号。 -q或–quiet或–silent 不显示任何信息。 -r或–recursive 此参数的效果和指定“-d recurse”参数相同。 -s或–no-messages 不显示错误信息。 -v或–revert-match 反转查找。 -V或–version 显示版本信息。 -w或–word-regexp 只显示全字符合的列。 -x或–line-regexp 只显示全列符合的列。 -y 此参数的效果和指定“-i”参数相同。 –help 在线帮助。

pattern正则表达式主要参数： \： 忽略正则表达式中特殊字符的原有含义。 ^：匹配正则表达式的开始行。 $: 匹配正则表达式的结束行。 <：从匹配正则表达 式的行开始。 >：到匹配正则表达式的行结束。 [ ]：单个字符，如[A]即A符合要求 。 [ - ]：范围，如[A-Z]，即A、B、C一直到Z都符合要求 。 . ：所有的单个字符。 *： 所有字符，长度可以为0。

4.grep命令使用简单实例 $ grep ‘test’ d* 显示所有以d开头的文件中包含 test的行。 $ grep ‘test’ aa bb cc 显示在aa，bb，cc文件中匹配test的行。 $ grep ‘[a-z]{5}’ aa 显示所有包含每个字符串至少有5个连续小写字符的字符串的行。 $ grep ‘w(es)t.\1′ aa 如果west被匹配，则es就被存储到内存中，并标记为1，然后搜索任意个字符(.)，这些字符后面紧跟着 另外一个es(\1)，找到就显示该行。如果用egrep或grep -E，就不用”\”号进行转义，直接写成’w(es)t.*\1′就可以了。

5.grep命令使用复杂实例 假设您正在’/usr/src/Linux/Doc’目录下搜索带字符 串’magic’的文件： $ grep magic /usr/src/Linux/Doc/* sysrq.txt:* How do I enable the magic SysRQ key? sysrq.txt:* How do I use the magic SysRQ key? 其中文件’sysrp.txt’包含该字符串，讨论的是 SysRQ 的功能。 默认情况下，’grep’只搜索当前目录。如果 此目录下有许多子目录，’grep’会以如下形式列出： grep: sound: Is a directory 这可能会使’grep’ 的输出难于阅读。这里有两种解决的办法： 明确要求搜索子目录：grep -r 或忽略子目录：grep -d skip 如果有很多 输出时，您可以通过管道将其转到’less’上阅读： $ grep magic /usr/src/Linux/Documentation/* | less 这样，您就可以更方便地阅读。

有一点要注意，您必需提供一个文件过滤方式(搜索全部文件的话用 *)。如果您忘了，’grep’会一直等着，直到该程序被中断。如果您遇到了这样的情况，按 ，然后再试。

下面还有一些有意思的命令行参数： grep -i pattern files ：不区分大小写地搜索。默认情况区分大小写， grep -l pattern files ：只列出匹配的文件名， grep -L pattern files ：列出不匹配的文件名， grep -w pattern files ：只匹配整个单词，而不是字符串的一部分(如匹配’magic’，而不是’magical’)， grep -C number pattern files ：匹配的上下文分别显示[number]行， grep pattern1 | pattern2 files ：显示匹配 pattern1 或 pattern2 的行， grep pattern1 files | grep pattern2 ：显示既匹配 pattern1 又匹配 pattern2 的行。

grep -n pattern files 即可显示行号信息

grep -c pattern files 即可查找总行数

这里还有些用于搜索的特殊符号： < 和 > 分别标注单词的开始与结尾。 例如： grep man * 会匹配 ‘Batman’、’manic’、’man’等， grep ‘<man’ * 匹配’manic’和’man’，但不是’Batman’， grep ‘<man>’ 只匹配’man’，而不是’Batman’或’manic’等其他的字符串。 ‘^’：指匹配的字符串在行首， ‘$’：指匹配的字符串在行 尾， 突然感到心累，grep命令好复杂T_T

今天就到这吧，后背酸痛，明天继续0.0