ANNOROAD0922 생물정보학 연습문제1, 데이터/newBGISeq5001.fq 및 data/newBGISeq5002. fq에서 BGISeq500 측정 플랫폼을 바탕으로 하는 진핵생물 유전자 그룹 DNA의 PE101 측정 데이터로 삽입 부분의 길이는 450bp이다.이 게놈 크기는 약 6M 정도로 알려져 있다.1) 이 테스트의 PE reads 수를 몇 쌍으로 집계해 보십시오.이론적으로 유전자 그룹의 99퍼센트 이상의 구역을 적어도 40X에 이르게 할 수 있습니까?추리와 계산의 과정과 결과를 간단하게 쓰십시오. 수치 계산은 R 등의 도구를 사용할 때 사용하는 코드를 쓰십시오.코드: 1, wc -l data/newBGIseq500_1.fq|awk '{print $1/4}' 159999 2, 참조 1

이전 단계의 결과에 근거하여 모든 데이터량을 계산하다

base

평균 깊이를 계산하면 깊이가 파송 분포(이론상)에 부합되고 평균 깊이는 그 기대에 부합된다.

dep

현재 상황을 계산합니다. 깊이가 40X보다 큰 구간의 백분율을 계산합니다.

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

[1] 0.9650074
이론적으로 유전자 그룹의 99퍼센트 이상의 구역을 적어도 40X에 이르게 할 수 없다
참조 2:

genome 

참조 3:

이전 단계의 결과에 근거하여 모든 데이터량을 계산하다

base

평균 깊이를 계산하고 (이론적으로 (가장 간단하게 이해하고 복종하는 분포 상황을 고려하지 않는다) 평균 깊이는 바로 그의 기대이다.

dep40

이론적으로 유전자 그룹의 99퍼센트 이상의 구역을 적어도 40X에 이르게 할 수 있다
결제 사고방식: 1. fq는 네 줄마다 하나의reads 정보를 표시한다.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(){         
next if(/^$/||/^>/);        
 	my $line=$_;          
 	my ($name,@seq)=split /
/,$line;         
 	my $list=join("",@seq);        
  	my $seqname=(split/\s+/,$name)[0];        
  	 my $length=length($list);        
   	 $hash{$seqname}=$length; 
    #       print "$seqname\t$length
"; 
    } 
    foreach my $key (sort { $hash{$b} <=> $hash{$a} } keys %hash ){         
    	print "$key\t$hash{$key}
";  
    } 

실행

perl $0 >length.txt

길이 작도: 위perl에서 생성한length를 연결합니다.txt.

pdf("length.pdf") 
lens 

예 1:

> qual  pdf("Length.pdf") 
> qual qual  qual$per  aa  aa  aa$per  sum  accu  for(i in 1:nrow(aa)){ sum = sum + aa[i,2] ，  accu[[i]] = sum } 
> aa$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

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=$total/2)&&($half==0)){print "N50: $x[$j]
";$half=$x[$j]}elsif($count>=$total*0.9){print "N90: $x[$j]
";exit;}}' output.scafSeq
  N50: 40176 
  N90: 3782

2. 시험 참고 디렉터리에 있는 파일 데이터/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  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: 7661779-7687550 이 영역에 포함된 변위 지점의 스크립트:

#!/usr/bin/perl 
use strict; 
use warnings; 
 my @sample; 
 my %hash;
  open (IN,"chr17.vcf") or die $!; 
  while () { 	
  	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

퀴즈 요구 사항: 퀴즈 디렉터리에 'exam' 이라는 디렉터리를 새로 만들고 그 안에 '1 SOAPdenovo', '2 Trio' 두 개의 하위 디렉터리를 만들어서 이 문제의 해답을 순서대로 하위 디렉터리에 저장하십시오.필요한 소프트웨어: SOAPdenovo, Jellyfish, VCF 조작 도구, R 등은 스스로 설치하세요.문제 풀이 사고방식과 정보 조회, 참고 시퀀스, 소프트웨어 다운로드 주소 등 스크립트가 아닌 절차를 Result에서 실행하십시오.txt 파일에서 설명, 실행 가능한 명령은 01work.sh、02_work.sh에서 다른 프로그램, 스크립트, 출력 파일 등은 필요에 따라 명명됩니다.