欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

生物信息学练习题-亚磊

程序员文章站 2024-03-04 12:08:53
...

ANNOROAD0922
生物信息学练习题
一、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

参考3：

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

base<-1599999*200

计算平均深度，（理论上(最简单的理解，不考虑服从的分布情况)），平均深度即其的期望

dep<-base/（6000000*0.99）
dep >40

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

结题思路：
1、 fq每四行表示一条reads的信息，1.fq和2.fq是成对存在的
2、 通过PE101、reads对数计算得到总碱基数，与6M基因组大小的99%的区域40X以上需要的碱基数作比较即可

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

下载地址：
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

提取序列长度：
脚本内容：

#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

长度作图：
接上面perl生成的length.txt。

pdf("length.pdf") 
lens <- read.table("length.txt") 
pdata <- cumsum(lens[,2]) 
plot(pdata,type="l",ylab='Total length',xlab = 'Seq num') 
dev.off()

示例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()

3）计算组装结果的N50。

 perl -e 'my ($len,$total)=(0,0);my @x;while(<>){if(/^[\>\@]/){if($len>0){$total+=$len;[email protected],$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基因上是变异情况进行输出，并说明变异位点的数目以及各样品的情况（纯合、杂合位点数目）。

/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
该区域中包含的变异位点的脚本：

#!/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;		
   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中，其他程序、脚本、输出文件等按需要命名。

相关标签： 0922考试 N50 组装 泊松分布 测序深度

上一篇： Mac 上使用 TeX

下一篇： java计算开方

推荐阅读