生物資訊練習###
阿新 • • 發佈：2018-12-11
一、data/newBGIseq500_1.fq和data/newBGIseq500_2.fq中是基於BGIseq500測序平臺的一種真核生物基因組DNA的PE101測序資料，插入片段長度為450 bp；已知該基因組大小約在6M左右。
1）請統計本次測序的PE reads數是多少對reads？理論上能否使基因組99%以上的區域達到至少40X覆蓋？請簡要寫出推理和計算的過程與結果，數值計算使用R等工具時請寫出所用程式碼。
2）請下載並安裝SOAPdenovo軟體，設定-K引數為35對該資料進行de novo組裝，並畫出組裝結果序列從長到短的長度累積曲線圖；
3）計算組裝結果的N50。
（1）perl指令碼
#!/usr/bin/perl -w use strict; my $file = "newBGIseq500_2.fq"; open(FH,$file) or die $!; my $reads = 0; while(<FH>){ if($_=~/^CL/){$reads+=1;} } close FH; print"reads數是：$reads\n"; my $depth=40; my $sum=6000000; my $coverage=0.99; my $length=200; my $final=$sum*$coverage*$depth/$length; print"一共需要reads數目為：$final\n"; if($reads>$final) {print"可以達到\n";} else{print"不能達到\n";}
（2）soapdenovo
1./home/stu8/BGI/software/SOAPdenovo2-master/SOAPdenovo-63mer all -s lib.cfg -K 35 -D 1 -o species>>ass.log lib.cfg max_rd_len=100 [LIB] avg_ins=450 reverse_seq=0 asm_flags=3 rank=1 pair_num_cutoff=3 map_len=32 q1=/home/stu8/BGI/exam/data/newBGIseq500_1.fq q2=/home/stu8/BGI/exam/data/newBGIseq500_2.fq 2.perl 1-2.pl species.scafSeq #usr/bin/perl -w use strict; $/=">"; my %hash; open(IN,"species.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"; } 3.Rscript 1-2.R #usr/bin/Rstript pdf("length.pdf") data<-read.table("qqq",head=F) colnames(data)<-c("name","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,type="l",xaxt="n",yaxt="n",main="length frequence",xlab="length",ylab="rate") axis(side=2,seq(from=0,to=1,by=0.1)) axis(side=1,seq(from=0,to=300000,by=50000),) grid(nx=25,ny=25) #length(mulData$length) #length(mulData$num/sum) dev.off()
(3)less species.scafStatistics
二、考試參考目錄下檔案data/chr17.vcf.gz，中是某trio家系的17號染色體的變異集合，參考序列為hg38。
1）編寫指令碼或選擇適當工具，統計vcf中變異位點的Qual值分佈情況，並畫圖展示。
2）選擇合適的工具或方法提取該家系在TP53基因上是變異情況進行輸出，並說明變異位點的數目以及各樣品的情況（純合、雜合位點數目）。
（1）
1.#統計 vcf 檔案中 Qual 值分別情況
less -N chr17.vcf | grep -w '^chr17' | awk -F "\t" '{print $2"\t"$6}' > chr17.qual
或者
vcftools --gzvcf chr17.vcf.gz --site-quality --out Q2.Rscript qual.R
#qual.R
pdf("QuallStatistics.pdf")
a <- read.table("chr17.qual",header = FALSE)
q <- a[,2]
hist(q,main = "Qual Hist")
dev.off()
（2）
1.
```/bin/vcftools --gzvcf chr17.vcf.gz --chr chr17 --from-bp 7668402 --to-bp 7687550 --het -c | head
or
zcat chr17.vcf.gz | grep '^chr17' | awk -F "\t" '{if ($2>=7668402 && $2<=7687550){print $0}}' | wc -l結果：
INDV O(HOM) E(HOM) N_SITES F
27DMBDM4YT 24 16.0 24 1.00000
7XKZJA3JWX 0 16.0 24 -2.00000
APRDKV0LDS 24 16.0 24 1.00000
2.
···/bin/vcftools --gzvcf chr17.vcf.gz --chr chr17 --from-bp 7668402 --to-bp 7687550 --recode --out tp533.sh 2-2.sh > 2-2.txt
for i in {10..12};do
cat tp53.recode.vcf | grep '#CHROM' |awk -F "\t" '{print $'$i'}'
echo -e "hom\t `cat tp53.recode.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1==$2){print $0}}' | wc -l`"
echo -e "het\t `cat tp53.recode.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1!=$2){print $0}}' | wc -l`"
done結果：
27DMBDM4YT
hom 9 #純合變異
het 8 #雜合變異
7XKZJA3JWX
hom 6
het 33
APRDKV0LDS
hom 8
het 8