生物信息学练习题

生物信息学练习题

一、data/newBGIseq500_1.fq和data/newBGIseq500_2.fq中是基于BGIseq500测序平台的一种真核生物基因组DNA的PE101测序数据,插入片段长度为450 bp;已知该基因组大小约在6M左右。

1) 请统计本次测序的PE reads数是多少对reads?理论上能否使基因组99%以上的区域达到至少40X覆盖?请简要写出推理和计算的过程与结果,数值计算使用R等工具时请写出所用代码。

代码:

1、 wc -l data/newBGIseq500_1.fq|awk ‘{print $1/4}’

1599999

2、

参考1

# 根据上一步结果计算全部的数据量

base<-1599999*200

# 计算平均深度,深度符合泊松分布(理论上),平均深度即其的期望

dep<-base/6000000

# 计算当前情况,深度大于40X的区间的百分比。

print (ppois(40,lambda=dep,lower=FALSE))

[1] 0.9650074

理论上不能使基因组99%以上的区域达到至少40X覆盖

参考2:

genome <- 6e6

read_length = 100

number_reads <- 6399996/2

depth = read_length * number_reads / genome

print(1- ppois(39 , depth))

[1] 0.975145

结题思路:

1、 fq每四行表示一条reads的信息,1.fq和2.fq是成对存在的

2、 通过PE101、reads对数计算得到总碱基数,与6M基因组大小的99%的区域40X以上需要的碱基数作比较即可

2) 请下载并安装SOAPdenovo软件,设置-K参数为35对该数据进行de novo组装,并画出组装结果序列从长到短的长度累积曲线图;

clip_image002

下载地址:

https://sourceforge.net/projects/soapdenovo2/files/latest/download

安装:

make

配置文件:

#maximal read length

max_rd_len=100

[LIB]

#average

insert size avg_ins=450

#if sequence needs to be reversed

reverse_seq=0

#in which part(s) the reads are used

asm_flags=3

#use only first 100 bps of each read

rd_len_cutoff=100

#in which order the reads are used while scaffolding

rank=1

# cutoff of pair number for a reliable connection (at least 3 for short insert size)

pair_num_cutoff=3

#minimum aligned length to contigs for a reliable read location (at least 32 for short insert size)

map_len=32

#a pair of fastq file, read 1 file should always be followed by read 2 file

q1=/home/stu27/data/newBGIseq500_1.fq

q2=/home/stu27/data/newBGIseq500_2.fq

命令执行

/home/stu27/soapdenovo/SOAPdenovo2-master/SOAPdenovo-63mer all -s /home/stu27/test/soapdenovo/example.config -K 35 -R -o output 1>ot1.log 2>ou1.err

提取序列长度:

脚本内容:

awk ‘!(NR%4)’ test.fq | awk ‘{print length($1);}’ > length

R

data=read.table(‘length’,header=FALSE)

data=as.matrix(data)

data=sort(data,decreasing=TRUE)

data=t(t(data))

qsum={}

apply(data,1,function(x)rbind(qsum,sum(data[1:x])))

pdf(“length.pdf”)

lines(qsum)

dev.off()

#usr/bin/perl -w

use strict;

$/=”>”;

my %hash;

open(IN,”output.scafSeq”) or die $!;

while(<IN>){

next if(/^$/||/^>/);

my $line=$_;

my ($name,@seq)=split /\n/,$line;

my $list=join(“”,@seq);

my $seqname=(split/\s+/,$name)[0];

my $length=length($list);

$hash{$seqname}=$length;

# print “$seqname\t$length\n”;

}

foreach my $key (sort { $hash{$b} <=> $hash{$a} } keys %hash ){

print “$key\t$hash{$key}\n”;

}

执行

perl $0 >length.txt

长度作图:

示例1:

> qual <- read.table(“length.list”)

> pdf(“Length.pdf”)

> qual<- ftable(qual)

> qual <- as.data.frame(qual)

> qual$per <- qual$Freq/sum(qual$Freq)

> aa <- sort(qual[,1],decreasing = T)

> aa <- as.data.frame(aa)

> aa$per <- as.data.frame(rev(qual[,3]))

> sum <- 0

> accu <- list()

> for(i in 1:nrow(aa)){ sum = sum + aa[i,2] , accu[[i]] = sum }

> aa$accu <- t(as.data.frame(accu))

> plot(x = aa[,1],y = aa[,3],type=”o”,col =”red”,xlab =”length”,ylab =”percentage”,main=”accuWarning message:ld”)

> dev.off()

示例2:

png(“length.png”)

data<-read.table(“zzz”,head=F)

colnames(data)<-c(“length”)

pdata<-table(data$length)

pframe<-data.frame(as.numeric(names(pdata)),as.numeric(pdata))

colnames(pframe)<-c(“length”,”num”)

rankData<-pframe[order(pframe[,1],decreasing=T),]

mulData<-data.frame()

sum=0

for(i in 1:nrow(rankData)){ sum=sum+rankData[i,2] row=cbind(rankData$length[i],sum) mulData=rbind(mulData,row) }

colnames(mulData)<-c(“length”,”num”) plot(mulData$length,mulData$num/sum,main=”序列长度累积曲线图”,xlab=”长度”,ylab=”累计比例”,type=”l”)

axis(side=2,seq(from=0,to=1,by=0.1),) #axis(side=1,seq(from=36000,to=300000,by=10000),) dev.off()

/export/home/stu1/export/stu1/data/export/stu1/data

3)计算组装结果的N50。

perl -e ‘my ($len,$total)=(0,0);my @x;while(<>){if(/^[\>\@]/){if($len>0){$total+=$len;push@x,$len;};$len=0;}else{s/\s//g;$len+=length($_);}}if ($len>0){$total+=$len;push @x,$len;}@x=sort{$b<=>$a}@x; my ($count,$half)=(0,0);for (my $j=0;$j<@x;$j++){$count+=$x[$j];if(($count>=$total/2)&&($half==0)){print “N50: $x[$j]\n”;$half=$x[$j]}elsif($count>=$total*0.9){print “N90: $x[$j]\n”;exit;}}’ output.scafSeq

N50: 40176

N90: 3782

二、考试参考目录下文件data/chr17.vcf.gz,中是某trio家系的17号染色体的变异集合,参考序列为hg38。

1) 编写脚本或选择适当工具,统计vcf中变异位点的Qual值分布情况,并画图展示。

提取qual值

/home/stu27/vcftools/vcftools_0.1.13/bin/vcftools –gzvcf chr17.vcf.gz –out Qual –site-quality

#会生成每个位点变异的质量值

作图:

> data <- read.table(“Qual.lqual”,header = TRUE)

> pdf(“Qual.pdf”)

> hist(data[,3],main = “Qual Hist”)

> dev.off()

2) 选择合适的工具或方法提取该家系在TP53基因上是变异情况进行输出,并说明变异位点的数目以及各样品的情况(纯合、杂合位点数目)。

Vcf共有10列 第6列是QUAL 变异质量

纯合和杂合的判断 vcf文件第10列 sample 0/1 0/0 两个数字是否相同 0是和ref一致 1是和ALT一致 1/2 表示杂合突变 有ALT1和ALT2两个基因型(有多个基因型的话就多于10列)

Vcftools可以输出vcf文件中的变异情况,由于需要起始和终止碱基位置,所以要去ucsc table broswer查看TP53基因对应的确定位置。

http://www.omicsclass.com/article/360

筛选vcf中指定位置的变异位点

/home/stu27/vcftools/vcftools_0.1.13/bin/vcftools –gzvcf chr17.vcf.gz –chr chr17 –to-bp 7687550 –out TP53 –recode –from-bp 7661779

以上语句输出的是包含变异位点的行

TP53位置信息:

Chromosome 17: 7,661,779-7,687,550

#统计变异数目和纯合杂合

变异位点数目:vcf文件 ref-alte 如果不为0 为变异位点 输出

纯合杂合:

取第10,11,12列

每列

先用冒号分割 取第一部分 再用/分隔 比较 1 和 2 列是否相同 相同纯合数加1 不同杂合数+1

该区域中包含的变异位点的脚本:

#!/usr/bin/perl

use strict;

use warnings;

my @sample;

my %hash;

open (IN,”chr17.vcf”) or die $!;

while (<IN>) {

chomp; #删除默认的分隔符

next if (/^##/);

my $line=$_;

if ($line=~/^#CHROM/) { (undef,undef,undef,undef,undef,undef,undef,undef,undef,@sample)=split/\t/,$line; #将每行数据以\t分隔,分别赋给变量

next;

}

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 27DMBDM4YT 7XKZJA3JWX APRDKV0LDS

my ($chrom,$pos,$id,$ref,$alt,$qual,$filter,$info,$format,@Other)=split/\t/,$line;

if ($pos>=7661779 && $pos <=7687550){

for (my $i=0;$i<@sample ;$i++) {

my $lis=(split/:/,$Other[$i])[0];

my ($aa,$bb)=split/\//,$lis;

if ($aa == $bb) {

if (exists $hash{$sample[$i]}{‘aa’}) {

$hash{$sample[$i]}{‘aa’} +=1;

}else{

$hash{$sample[$i]}{‘aa’}=1;

}

}else{

if (exists $hash{$sample[$i]}{‘ab’}) {

$hash{$sample[$i]}{‘ab’} +=1;

}else{

$hash{$sample[$i]}{‘ab’}=1;

}

}

}

}

}

close IN;

print “#Sample\tHomozygous\tHeterozygote\n”;

foreach my $name (keys %hash) {

print “$name\t$hash{$name}{‘aa’}\t$hash{$name}{‘ab’}\n”;

}

#Sample Homozygous Heterozygote

27DMBDM4YT 53 12

APRDKV0LDS 53 12

7XKZJA3JWX 13 52

答题要求:

请在答题目录下新建一个名为“exam”的目录,并在其中建立2个子目录: “1_SOAPdenovo”、“2_Trio”,将该题目的解答依次保存在子目录中。

需要用到的软件:SOAPdenovo、Jellyfish、VCF操作工具、R等等请自行安装。

请将解题思路与信息查询、参考序列、软件下载地址等非脚本运行的步骤在result.txt文件中进行说明,可运行命令保存在01_work.sh、02_work.sh中,其他程序、脚本、输出文件等按需要命名。