搜索
查看: 2529|回复: 0

[Other] Biostar: 课程9、10-测序文件的质控&在线序列比对

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-4-13 19:45:53 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:https://mp.weixin.qq.com/s/OVd-2T_6_4EgmjPe3oxZCg


翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

#为作者的注释
#*为译者的注释

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

# 到正则表达式测试网站去实战一下。
# 在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两两比对的网站来进行比对。

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

>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:[url=]http://www.cnblogs.com/tsw1107/p/2264a01aeec481d2044dfeda01417c64.html[/url]

欢迎到微信公众号订阅我们
生信媛
bio_sxy




上一篇:求助:Delly 做Somatic SV calling 的时候数据输出不完整。
下一篇:TCGA中miRNA数据的提取
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-19 01:19 , Processed in 0.044738 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.