|
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./project/bowtie-bio/bowtie2/2.2.4/bowtie2-2.2.4-macos-x86_64.zip
# Linux上用： #curl -OL http://downloads./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.fqREAD2=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.bamsamtools sort -f tmp.bam $name.bam
done
# 删除中间文件。
rm -f tmp.sam tmp.bam
# 修改名字，使得他们不再含有两个后缀。
mv bwa.sam.bam bwa.bammv bow.sam.bam bow.bam
# 对数据进行索引。
for name in *.bam
do
samtools index $name
done
# 删除这个程序生成的所有文件。
# rm -f *.bam *.bai *fq *.txt *.sam *.gff
参考资料1：https://www./article/655.html
|
