Biostar：课程17、18

查看原文

其他

Biostar：课程17、18

原创 2017-08-10 István Albert 生信媛

翻译：生物女博士

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

#为作者的注释

#*为译者的注释

Lecture 17 - 用awk进行简单的编程

# 用awk进行简单的编程

# 我们之前已经做过这个了，但是为了让大家我们了解我们正在做什么，我们来再做一遍。

# 提取埃博拉基因组的编码特征、基因以及编码序列。

efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb
readseq -format=GFF -o NC.gff NC.gb

#*可以跳过上边两步，使用拟南芥gff文件。

#*curl -O ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

#*我们可以先来查看一下gff格式是什么样子的

less -S TAIR10_GFF3_genes.gff |head

#*gff文件是tab分隔的文件。第1列是染色体信息；第2列是gff注释数据来源；第3列为特征（feature）即属于gene还是mRNA还是CDS等等；第4和5列分别是这个特征序列的起始和终止位置，第6列是得分，可以是序列相似性比对时的E-values值或者基因预测是的P-values值。”.”表示为空（来自参考资料1）；第7列是表示序列的方向：正义链为+，反义链为-；第8列仅为对CDS的注释，表示起始编码的位置，有效值为0、1、2（来自参考资料1）；第9列为注释信息。

# 找到每个特征的长度，注意“魔法”。

#*如果使用的是拟南芥的gff文件为操作对象，请自行将之后所有命令的NC.gff改为TAIR10_GFF3_genes.gff

#*$1, $2, $3指的就是第1-3列的数据

cat NC.gff | awk ' { print $1, $2, $3 } ' | head -5

# 这（基本上）等同于截取列。

cat NC.gff | cut -f 1,2,3 | head -5

# 我们可以进行计算。每个特征序列长度是多少？

#*第5列的数值减去第4列的数值后+1，即得到特征序列的长度

cat NC.gff | awk ' { print $3, $5-$4 + 1 } ' | head -5

# 我们可以利用模式匹配来提取CDS特征。

cat NC.gff | awk '$3 =="gene" { print $3, $5-$4 + 1, $9 } '

# 计算所有基因的累积长度。

cat NC.gff | awk '$3 =="gene" { len=$5-$4 + 1; size += len; print "Size:", size } '

# 计算所有CDS的累积长度。

cat NC.gff | awk '$3 =="CDS" { len=$5-$4 + 1; size += len; print "Size:", size } '

# 基因组有多大？有很多方法可以知道，比如从SAM文件的头文件（header）。

samtools view -h ../lec16/results.bam | head -2

# 基因组上有多少是基因。

#*如果是拟南芥，把perc=size/18959对应改成perc=size/119667750，119667750是拟南芥（Col-0）基因组的大小。

#*可以用下边的代码自行计算。

#*cat TAIR10_GFF3_genes.gff |awk '$3 == "chromosome"{len=$5-$4 + 1; size += len; print "Size:", size } '

cat NC.gff | awk '$3 =="CDS" { len=$5-$4 + 1; size += len; perc=size/18959; print "Size:", size, perc } '

# 我们需要使awk只用tab分隔符进行分隔。

# 用空格符分隔虽然有用，但是常常导致烦人的bugs

# 告诉awk，以tab作为F（字段分隔符）和OFS（输出字段分隔符）

# 把下边这行命令放到 .profile 或 .bashrc 文件中

alias awk="awk -F '\t' -v OFS='\t'"

# 使设置生效。

source ~/.profile

# 生成新的含有基因和CDS的gff文件。

# The first line specifies what kind of file this is.

head -1 NC.gff > NC-genes.gff
head -1 NC.gff > NC-cds.gff

# 根据特征（features）把文件分开。

cat NC.gff | awk ' $3=="gene" { print $0 }' >> NC-genes.gff
cat NC.gff | awk ' $3=="CDS" { print $0 }' >> NC-cds.gff

# Sam文件是tab分隔符的，可以用awk来处理。

# 上周的数据中，有多少碱基被覆盖了超过50次？

samtools depth ../lec16/results.bam | awk '$3 > 50 { print $0 } ' | wc -l

# 有多少的模板长度超过50 bp。

samtools view ../lec16/results.bam | awk ' $9 > 50 { print $0 } ' | wc -l

从GFF文件中分离提取基因名字。

$3 == "gene" {
# 通过 ; 分离提取第9列

split($9, x, ";")
# 基因名字是第一个元素。

# 通过空格符切分。基因名字是第二个元素

split(x[1], y, " ")

# 去除基因名字旁边的双引号

name = y[2]

#*对于拟南芥的gff，通过“=”切分。基因名字是第二个元素：
#*split(x[1], y, "=")
#*name = y[2]# 用空格全局替换 " 符号

# 由于 " 是一个特殊字符，我们必须写成 \"

#*反斜杠\表示转义符。

gsub("\"", "", name)

# 打印特征类型、基因名字以及大小。

print $3, name, $5 - $4 + 1}

#*最后，我们可以写成下边这条命令

#*cat TAIR10_GFF3_genes.gff |awk '$3 == "gene"{split($9,x,";");split(x[1], y, " ");name = y[2];gsub("\"", "", name); print $3, name, $5 - $4 + 1 } '

Lecture 18 - 序列比对工具的对比

# 下载并安装比对软件bowtie2

cd ~/src

# Mac OSX上用：

curl -OL http://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.2.4/bowtie2-2.2.4-macos-x86_64.zip

# Linux上用：

#curl -OL http://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.2.4/bowtie2-2.2.4-linux-x86_64.zip

# 这个已经是个二进制文件了，所以可以直接执行。

unzip bowtie2-2.2.4-macos-x86_64.zip

# 把比对软件以及相关程序链接到bin文件夹。

ln -s ~/src/bowtie2-2.2.4/bowtie2 ~/bin/
ln -s ~/src/bowtie2-2.2.4/bowtie2-align ~/bin/
ln -s ~/src/bowtie2-2.2.4/bowtie2-build ~/bin/

# 对参考基因组建立索引。

# 这是bowtie2的索引。

bowtie2-build ~/refs/852/NC.fa NC.fa

# 把突变弄成一定的格式，使得可以在浏览器上展示。

# 把它变为gff文件。

#*mutations.txt这个文件是在lecture15中生成的。并且存放在目录~/edu/lec15下。

cat mutations.txt | awk ' {print $1, "wgsim", "mutation", $2, $2, ".", "+", ".", "." }' > mutations.gff

# 运行比对程序。

#*脚本见文末

bash compare.sh

# 修改脚本使之生成有较大的错误率的fastq文件。

# 当你修改了错误率后，计算一下比对上的reads数目。

samtools view -cF 4 bwa.bam
samtools view -cF 4 bow.bam

最终脚本compare.sh

# 对比两个比对软件的输出。
# 使用方法: bash compare.sh
# 在出错的地方停止运行。
set -ue

# 数据文件名。
READ1=r1.fq
READ2=r2.fq

# 参考序列文件。两个比对软件需要分别对它进行索引。
REFS=~/refs/852/NC.fa

# 从基因组中生成模拟reads。
# 编辑错误率并重新运行这个pipeline。
wgsim -N 10000 -e 0.1 $REFS $READ1 $READ2 > mutations.txt

# 将mutations文件转成GFF文件。
cat mutations.txt | awk -F '\t' ' BEGIN { OFS="\t"; print "##gff-version 2" } { print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "name " $3 "/" $4 }' > mutations.gff

# 我们需要添加一个read组到mapping中。

GROUP='@RG\tID:123\tSM:liver\tLB:library1'

# 运行bwa并创建bwa.sam文件。
bwa mem -R $GROUP $REFS $READ1 $READ2 > bwa.sam

# 运行bowtie2并创建bow.sam文件。
bowtie2 -x $REFS -1 $READ1 -2 $READ2 > bow.sam

# 调整bowtie2# bowtie2 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x $REFS -1 $READ1 -2 $READ2 > bow.sam
# 对每个sam文件，都转换成bam（sam的二进制）格式。
for name in *.sam;
do
samtools view -Sb $name > tmp.bam
samtools sort -f tmp.bam $name.bam
done

# 删除中间文件。
rm -f tmp.sam tmp.bam

# 修改名字，使得他们不再含有两个后缀。
mv bwa.sam.bam bwa.bam
mv bow.sam.bam bow.bam

# 对数据进行索引。
for name in *.bam
do
samtools index $name
done

# 删除这个程序生成的所有文件。
# rm -f *.bam *.bai *fq *.txt *.sam *.gff