生物信息学练习题-赞哥

ANNOROAD0922

一、data/newBGIseq500_1.fq 和 data/newBGIseq500_2.fq 中是基于 BGIseq500 测序平台的

一种真核生物基因组 DNA 的 PE101 测序数据，插入片段长度为 450 bp；已知该基因组大小
约在 6M 左右。

1）请统计本次测序的 PE reads 数是多少对 reads？理论上能否使基因组 99%以上的区域达到

至少 40X 覆盖？请简要写出推理和计算的过程与结果，数值计算使用 R 等工具时请写出所用代码。

wc -l newBGIseq500_2.fq| awk '{print($1/4)}'

1599999

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

#推导过程如下
1）先计算平均深度，平均深度为53.3X，R代码如下：

base <- 1599999*200
dep <- base/6000000
dep
[1] 53.3333

2）根据正态分布的68-95-99 法则
距离均值3个标准差的区域围起来的范围占到了总体的99%
Hiseq平台的重测序的经验 标准差为10X
要使基因组99%的区域至少达到40X，需要的最低平均深度= 40X + 10*3 X
既平均深度至少要达到70X才行，本次测序平均深度为53.3X，所以答案为 不可能

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

SOAPdenovo-63mer all -s example.config -K 35 -R -o output 1>ot1.log 2>ou1.err

scafold_lenth.pl

#usr/bin/perl -w
use strict;
$/ = ">";
open(IN,"output.scafSeq") || die $!;
while (<IN>){
 	chomp;
 	next if (/^$/);
 	my @tmp = split /\n/,$_,2;
 	$tmp[1] =~ s/\n//g;
 	print length($tmp[1])."\n"
}
close IN;

perl scafold_lenth.pl |sort -k 1rn |awk 'BEGIN{sum=0}{if(!NF){next};sum+=$1;print sum}' >len.txt

R 画图

lens <- read.table("len.txt")
plot(x=seq(1,length(lens$V1)),y=lens$V1,type="l")

3）计算组装结果的 N50。

先通过上一步的到的len.txt看最后一行得到scaffold的总长度 为 6181532
再用R计算总长度的一半

> 6181532/2
[1] 3090766

cat len.txt |awk '{if($1>=3090766){print $1;exit}}'

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

考序列为 hg38。

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

awk '{if(/^#/){next};print $6}' chr17.vcf >qual.txt

R实现画图

pdf("qual.pdf")
qual <- read.table('qual.txt')
hist(qual$V1,main = "Qual Hist")
dev.off()

2）选择合适的工具或方法提取该家系在 TP53 基因上是变异情况进行输出，说明变异位

点的数目以及各样品的情况（纯合、杂合位点数目）。
通过ensembl 查到位置信息： 7,661,779-7,687,538

vcftools --vcf chr17.vcf --chr chr17 --to-bp 7687538  --out TP53 --recode --from-bp 7661779
grep -v "^#" TP53.recode.vcf |wc -l

grep -v "^#" TP53.recode.vcf |cut -f 10|awk -F ":" '{print $1}'|awk -F "/" 'BEGIN{HO=0;HE=0}{if ($1==$2){HO+=1}else(HE++)}END{print "HO:"HO,"HE:"HE}’

grep -v "^#" TP53.recode.vcf |cut -f 11|awk -F ":" '{print $1}'|awk -F "/" 'BEGIN{HO=0;HE=0}{if ($1==$2){HO+=1}else(HE++)}END{print "HO:"HO,"HE:"HE}'

grep -v "^#" TP53.recode.vcf |cut -f 12|awk -F ":" '{print $1}'|awk -F "/" 'BEGIN{HO=0;HE=0}{if ($1==$2){HO+=1}else(HE++)}END{print "HO:"HO,"HE:"HE}'