Biostar课程：5、6

#为作者的注释

#*为译者的注释

Lecture 5 - 编写可重复使用的脚本。获取子序列。

# 通过locus号来获取一条序列。

efetch -db nucleotide -id KM233090 -format fasta > KM233090.fa

# 这条序列也可以通过accession号来获取。

efetch -db nucleotide -id 667853062 -format fasta > 667853062.fa

# 可以一次性获取多个id下的序列。

efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta > multi.fa

# 查看文件有多少行。

wc -l multi.fa

# 数有多少个头文件。

cat multi.fa | grep ">" | wc -l

# 查看头文件。

cat multi.fa | grep ">"

# 我们也可以从数据中只获取一个子集。

efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta -seq_start 1 -seq_stop 10 > multi.fa

# 我们开始来获取所有的埃博拉病毒基因组的开头。

#*这条命令获取的是每条序列的第1-10个碱基。

esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10 > starts.fa

# 每次都手动输入命令太费事，是时候使点黑技能来解决这个问题了。开始写shell脚本吧。

# shell脚本就是很多命令的集合，是可以重复运行的。

# 把下边的这些行放到一个文件中，通常以.sh为后缀。我将其命名为get_seq.sh

# esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10

#现在，你可以把一个脚本里的所有命令运行了。

bash get_seq.sh > starts.fa

# 能从外部指定这个脚本的某个部分，会更好。

# 两个$符号将被来自命令的参数取代。 详情参见下边的内含文件。

bash get_seq.sh 1 10 > starts.fa

# 现在，只把这个文件中的序列提取出来。

# -A 1 是指把匹配行及其下一行打印出来。

#  -v 是指把不匹配的行打印出来。

# 不要直接一下子输入全部语句，而要每输入一个语句，检查结果，再输入第二个，如此反复。

cat starts.fa | grep ">" -A 1 | grep -v ">" | grep -v "-"

我们第一版的脚本：

#!/bin/sh## 从所有的序列中把核苷酸子序列抓取出来。# 存放在PRJNA257197 这个项目（project）里。## 预期有start和stop这两个参数。# 这个设置使脚本停止在错误和未定义变量上。set -ue# 通过Entrez Direct访问NCBI。esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start $1 -seq_stop $2

我们前面脚本的升级版。当你加上函数时，生成不同的命令就变得更为容易了：

# 从所有的序列中把核苷酸子序列抓取出来。# 存放在PRJNA257197 这个项目（project）里。# 预期有start和stop这两个参数。# 这个设置使脚本停止在错误和未定义变量上。set -ue

function esearch_ebola() {    # 搜索ebola（埃博拉）的核苷酸序列。
    esearch -db nucleotide -query PRJNA257197
}

function efetch_limits() {    # 你可以把变量表示什么变得更明确。
    # 对于更为复杂的命令，这应该会帮助。
    START=$1
    END=$2
    # 运行efetch.
    efetch -db nucleotide -format fasta -seq_start $START -seq_stop $END}# 把两个命令链起来获取数据。# 请留意着两个参数是如何传递到函数中的。esearch_ebola | efetch_limits $1 $2

Lecture 6 - 编码和fastq格式

# 使用python后跟着-c这个命令来运行字符串。

python -c 'print "Hello World"'

# 字符的序数（整数）。

python -c 'print ord("A")'

#  整数的字符(字母代码) 。

python -c "print chr(65)"

# 找到字符对应的Sanger（+33）编码质量值

python -c 'print ord("I") - 33'

# 找到字符对应的Illumina 1.5 (+64)编码质量值

python -c 'print ord("I") - 64'

# 找到40这个值对应的Sanger (+33)编码质量字符

python -c 'print chr(40 + 33)'

# 找到40这个值对应的Illumina 1.5 (+64)编码质量字符

python -c 'print chr(40 + 64)'

#* 关于这个质量值和字符的对应关系，在《小白生信学习记3》有详细介绍。

# 从这个课程的网址中获取一个测试数据集。

curl -O http://www.personal.psu.edu/iua1/courses/files/2014/illumina-data.fq.gz

#* 原文是curl -O http://localhost:8080/courses/files/2014/illumina-data.fq.gz

#* 这个地址是当时上课的学生用的。我们是不能这么用的。

# 你可以解压这个gzip文件 illumina-data.fq.gz ，也可以不解压，并用gzcat (Mac OSX)或者zcat (Linux)来打开。当你有很多的数据集时，牺牲读取速度来节省存储也是值得的。

time gzcat illumina-data.fq.gz > tmp

# 解压数据，同时保留原压缩文件。

# 没有 -c 的话，你可以无需重定向一个输出，它会在当前文件夹下解压。

gunzip -c illumina-data.fq.gz > illumina-data.fq
time cat illumina-data.fq > tmp

# 查看FASTQ记录 (Mac OSX)

gzcat illumina-data.fq.gz | head -4

# 查看FASTQ记录  (Linux)

zcat illumina-data.fq.gz | head -4

# 拆包一个SRA数据和机器下来的原始数据集

# 此处会生成2个文件。

fastq-dump --split-files SRR1553605

#*这个数据之前是下载过的。详情参照上两节的翻译。

# 显示其中一个fastq的记录。

cat SRR1553605_1.fastq | head -4

# 简单粗暴地检查一下质量，在前10,000行里是否有连续的###字符。

gzcat illumina-data.fq.gz | head -10000 | grep '###' |  wc -lcat SRR1553605_1.fastq | head -10000 | grep '###' | wc -l

# 哪个看起来要好些？

# 只看前面/最后的几行的质量是有些冒险的，因为数据的排序可能不是随机的。通常我们会生成一个随机样本（我们在后边会看到如何去实现）。

# 我们可以在数据里搜索pattern并着色。

# 这些都是起始密码子吗？是否存在pattern TATA?

gzcat illumina-data.fq.gz | head -10000 | grep 'ATG' --color=always | head

gzcat illumina-data.fq.gz | head -10000 | grep 'TATA' --color=always | head

# 有看起来像Illumina接头: GATCGGA的序列吗？

gzcat illumina-data.fq.gz | head -10000 | grep 'GATCGGA' --color=always | head

翻译：生物女博士

注：本系列课程翻译已获得原作者授权。