1. 实用技巧

1.1 日志与后台输出

  当我们的指令会输出一系列长内容时,我们会考虑用日志文件记录输出.

kingfisher get -r SRR6451531 -m aws-http aws-cp ena-ftp 1>download.log 2>download.err

  有时候,程序开发的时候会将部分标准输出也归类到错误输出里去,所以我们很少区分正常输出的$log$和$err$,而是直接用一个文件记录.

kingfisher get -r SRR6451531 -m aws-http aws-cp ena-ftp 1>download.log 2>&1

  但是更多的,我们是将它放在后台输出.如果电脑与服务器远程连接终止后,后台运行的程序会终止.

kingfisher get -r SRR6451531 -m aws-http aws-cp ena-ftp 1>download.log 2>&1 &

  如果想终止连接也不停止运行,这时候我们使用$nohup$.

nohup 英文全称 no hang up（不挂起），用于在系统后台不挂断地运行命令，退出终端不会影响程序的运行. nohup 命令，在默认情况下（非重定向时），会输出一个名叫 nohup.out 的文件到当前目录下，如果当前目录的 nohup.out 文件不可写，输出重定向到 $HOME/nohup.out 文件中.

nohup kingfisher get -r SRR6451531 -m aws-http aws-cp ena-ftp 1>download.log 2>&1 &

1.2 监控后台程序

  我们使用htop来监控后台进程.

htop -u ug1700

1.3 杀死进程

  可以看下面三种指令.

2. 批量下载数据

  如果我们只是下载一个数据这样写还比较轻松,但如果是多个数据,这样写无疑自寻死路.这里我们要去NCBI数据库下载meta数据.

  下载好后,将后缀改成可以用$excel$打开的形式,比如$tsv$我们将$sample$对应名与ID存入新建文件sra_info.tsv.

  我们用awk指令来生成批量下载数据的命令并传入文件run_kf.sh.

awk '{print "kingfihser get -r " $1 " -m aws-http aws-cp ena-ftp >"$1".log 2>&1 &"}' sra_info.tsv  > run_kf.sh

  我们将SRR6451531的两个fastq文件下载好.

3. 过滤和质控

  数据中存在不合格的数据,我们需要将此数据扔掉,这就是过滤.质控就是看看数据的质量(比如$Q30$的值是多少).严格来说,数据到手后应该经过过滤质控再过滤的环节.

3.1 如何过滤低质量的reads

  现在二代测序成本降低,我们过滤低质量的$reads$的方法就是将这条$reads$丢弃.那么丢弃的标准是什么呢?我们之前讲过碱基质量值,下面来自他人博客引用.

碱基质量中的质量非物理学中的“质量”，碱基质量值（quality score，Q-score），在生物物理学中是碱基识别出错概率的整数映射,公式是$Q=-10 \times lgP$,其中P为碱基识别出错的概率.简单点来说是，在下机数据中每个序列的每个碱基都有一个质量值信息，我们通过识别这个质量值信息就可以了解到这个碱基被识别出错的概率是多少。

  我们过滤低质量reads,有多个不同的角度.

  使用的是下图两个标准.低质量碱基阈值是指如果碱基质量低于此阈值,就认为是低质量碱基.低质量碱基比例阈值是指当低质量碱基的比例超过此数值,就认为该read是低质量read.

  当reads过短时,我们也认为是低质量reads.

  还有一个角度是$adapter$.首先我们需要阐述什么是$adapters$污染.$adapter$污染指的是由插入片段长度不够，测序仪读到的测序引物等序列.我们要把测到接头序列的reads全部扔掉.去除有两种方法,如下图.

• 指定序列.指定接头序列的具体碱基排列,然后$filter$会根据序列在reads里匹配,匹配上的就去掉.
• 基于$overlap$自动检测.我们测序得到的$reads$存在大量读到接头的情况.我们不断读到接头,这种大量相同的序列不断出现,很有可能就是$adapter$.

3.2 使用软件过滤

  我们用于过滤的包是$fastp$.安装后,根据官方文档使用.
  过滤的代码如下:

nohup fastp -i ../12.raw_data/SRR6451531_1.fastq.gz -I ../12.raw_data/SRR6451531_2.fastq.gz -o CS1_1_1.fq.gz -O CS1_1_2.fq.gz -h CS1_1.html -j CS1_1.json &

  在说明原理之前,我想先后知后觉地记录一下$fastq.gz$文件是什么.我们用$zcat$查看压缩文件时,直接看到的是fastq文件的内容.也就是说我们的fastq文件是直接压缩成gz文件的.我们将$fastq.gz$解压,获得的就是$fastq$文件.

  原理在官网已经阐述.这和我们之前阐述的3个角度一致.质量、比例、长度.但是要说明的是,官网也可以通过复杂度过滤.复杂度定义为与其下一个基数不同的基数的百分比$(base[i] != base[i+1])$,一般处于默认状态.

  我们过滤玩完后会得到两个文件,一个是$html$格式,一个是$json$格式,这两个主要是记录此时过滤的结果情况.见下图.

  关于$Gerernal$的一些项,还是说明一下.比如duplication rate,如果两条序列测出来的是一模一样的,那么就会判断为$duplicate$,这个在PCR里会比较高.在一般的DNA测序中,大概是$4$ ~ $8$%左右.在RNA测序会比较高,20%~30%.这是因为每个基因的表达量不同,如果相同基因表达量较高,被抽到的几率就会比较高.

  下图是过滤前后对比.双端测序将一次各自双端测序算作1个reads,因此这里total reads乘上150就算作数据量(一条被打碎的$DNA$片段 150碱基).
这里我不是很确定了,感觉要复习一下RNA测序,一个迷惑的点是老师say这里要乘300才是数据量,但实际是乘150

  感觉也不用纠结那么多,这里有个类似的提问.

  这里有RNA测序说的比较详细的,当时建库是只建了RNA方向的cDNA库Click.

  当我们要对不同样本批量生成指令时,也很简单.

awk '{print "nohup fastp -i ../12.raw_data/" $1 "_1.fastq.gz -I ../12.raw_data/" $1 "_2.fastq.gz -o " $2 "_1.fq.gz -O " $2 "_2.fq.gz -h " $2 ".html -j " $2 ".json &"}' ../12.raw_data/sra_info.tsv