Biostar: 课程9、10
翻译：生物女博士
注：本系列课程翻译已获得原作者授权。
#为作者的注释
#*为译者的注释
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：
欢迎关注我们