搜索
查看: 2198|回复: 0

[Other] Biostar课程:5、6-编写可重复使用的脚本

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

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


#为作者的注释

#*为译者的注释

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 ... 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




翻译:生物女博士

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

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







上一篇:IRanges说明(3/3)
下一篇:Biostar:7、8-二代测序质控:fastqc、prinseq、trimmomatic、seqtk
回复

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.