Biostar: 课程9、10

原创 István Albert 生信媛

翻译：生物女博士

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

#为作者的注释

#*为译者的注释

Lecture 9 - 质控：模式匹配和接头剪切

# 这个数据集是来自上一节课

fastq-dump --split-files SRR519926

# 注意：模式搜索发生在行上，如果是sequences wrap则需要别的方法（译者没有搞懂sequences wrap该怎么翻译，orz，如有好的翻译，求文末留言。）大多数仪器产生reads的方式是一行一条，所以这个方法适用。但记住，这也可能匹配到质量一行的字符串（虽然通常不太可能）。其中一个解决方案是先把序列提取出来。

# 搜索开头为ATG的行

#*在正则表达式里，行开头用^表示

cat SRR519926_1.fastq | egrep "^ATG" --color=always | head

#搜索行末为ATG的行

#*在正则表达式里，$表示匹配输入字符串的结束位置

cat SRR519926_1.fastq | egrep "ATG\$" --color=always | head

#*egrep表示扩展正则表达式，用法上grep –E等同于egrep（参考资料1）

#*上文中的"\"是转义的作用，这里去除或者保留并不影响匹配结果。

# 搜索TAATA或TATTA模式，这是一个字符范围

cat SRR519926_1.fastq | egrep "TA[A,T]TA" --color=always | head

# 搜索TAAATA 或者 TACCTA, 这些是一些分组

cat SRR519926_1.fastq | egrep "TA(AA,CC)TA" --color=always | head

#*译者测试了一下，上一段命令是无法匹配的。译者认为，分组时，不能","分开，应使用"|"

#*可以试一试把","改为"|"：

cat SRR519926_1.fastq | egrep "TA(AA|CC)TA" --color=always | head

#利用元字符批量匹配

# 搜索 “TA后有0个或多个A，接着是TA”

cat SRR519926_1.fastq | egrep "TA(A*)TA" --color=always | head

# 搜索“TA后有1个或多个A，接着是TA”

cat SRR519926_1.fastq | egrep "TA(A+)TA" --color=always | head

# 搜索“ TA后有2个或5个A，然后是TA”

cat SRR519926_1.fastq | egrep "TAA{2,5}TA" --color=always | head

# 到正则表达式测试网站去实战一下。

# http://regexpal.com/

# 在reads末端匹配IIlumina接头

# 在任何位置匹配AGATCGG，AGATCGG后可以跟着任何数目的碱基。

cat SRR519926_1.fastq | egrep "AGATCGG.*" --color=always | head

# 查看FastQC 和Trimmomatic的污染序列和接头序列文件。

ls ~/src/FastQC/Configuration/

# 文件内容涵盖一些已知的接头序列。

more ~/src/FastQC/Configuration/adapter_list.txt

# 文件内容涵盖已知的污染序列。

more ~/src/FastQC/Configuration/contaminant_list.txt

# 进行接头切除时，我们需要找到含有接头序列的文件。

# 你可以自己创建自己的接头文件，或者使用Trimmomatic自带的

# 我们来创建一个Illummina的接头文件吧。

echo ">adapter" > adapter.fa
echo "AGATCGGAAGAG" >> adapter.fa

# 我们来进行质量和接头的修剪。

trimmomatic SE SRR519926_1.fastq trimmed.fq ILLUMINACLIP:adapter.fa:2:30:10 SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:adapter.fa:2:30:5

# Trimmomatic自带了一些接头，我们可以来用一下。

# 依据文件内容不同，可能有不同的操作模式。

# 这是一个palindromic接头。 Trimmomatic通过接头名字来确定操作模式。

#*留意你的安装版本号。如果安装了更新版本，需要对应改一下路径名称。

ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-PE.fa

# 这是一个扩展的palindromic接头文件，它也列出了单个接头。

ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-PE-2.fa

# 这些是一些经典的接头。

ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-SE.fa

trimmomatic SE SRR519926_1.fastq trimmed.fq ILLUMINACLIP:adapter.fa:2:30:10 SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:TruSeq3-SE.fa:2:30:5

# 在双端（paired end）模式下运行时，两个测序文件需要被同时过滤和剪切。

# 命令行混乱随之而来。

trimmomatic PE SRR519926_1.fastq SRR519926_2.fastq trim_1P.fq trim_1U.fq trim_2P.fq trim_2U.fq ILLUMINACLIP:TruSeq3-PE-2.fa:2:10:10 SLIDINGWINDOW:4:30 TRAILING:20

# 一个稍微简单的选择是，加一个 -baseout

Lecture 10 - 比对

# 这节课不需要用到命令行工具。

# 我们会通过EBI两两比对的网站来进行比对。

# http://www.ebi.ac.uk/Tools/psa/

# 下边两个蛋白序列是从Marketa Zvelebil写的Understanding Bioinformatics里来的一个有趣的例子。

# http://www.amazon.com/Understanding-Bioinformatics-Marketa-Zvelebil/dp/0815340249

>1
THISLINE

# 第二个蛋白。

>2
ISALIGNED

# 你需要设定不同参数，来进行全局比对和局部比对。

# 如果有时间，我们来研究一下埃博拉病毒。

# 还记得那些病毒基因组的开头有时候看着差异很大么？

# 我们从两个差异最大的基因组中提取前100bp，并比对。

# KM034561 and KM233037

efetch -db nucleotide -id KM034561 -seq_start 1 -seq_stop 100 -format fasta > 1.fa

efetch -db nucleotide -id KM233037 -seq_start 1 -seq_stop 100 -format fasta > 2.fa

# 用Needle比对这两条序列。你发现什么了？

# 记住，CGAATAACT 是最常见的起始序列。

参考资料1：http://www.cnblogs.com/tsw1107/p/2264a01aeec481d2044dfeda01417c64.html

欢迎关注我们

继续滑动看下一个