生物資訊學練習題-亞磊
阿新 • • 發佈：2018-12-11
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組裝，並畫出組裝結果序列從長到短的長度累積曲線圖；
配置檔案：
#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中，其他程式、指令碼、輸出檔案等按需要命名。