周六考试（四）

00.data

cp /home/stu33/exam/1_SOAPdenovo/*.gz ./

01.filter

../software/SOAPnuke-master/SOAPnuke  filter -1 ../00.data/newBGIseq500_1.fq.gz -2 ../00.data/newBGIseq500_2.fq.gz -C newBGIseq500_clean_1.fq.gz -D newBGIseq500_clean_2.fq.gz -l 10 -q 0.1 -n 0.01 -G 1 -Q 2 -o ./ &

02.SOAPdenovo

cp ../software/SOAPdenovo2-master/example.config ./my.config

../software/SOAPdenovo2-master/SOAPdenovo-63mer all -s my.config -K 35 -o newBGIseq500 1>soap.log 2>soap.info

03.callSNP

../../software/bwa-0.7.17/bwa index -a bwtsw newBGIseq500.scafSeq
../../software/bwa-0.7.17/bwa mem newBGIseq500.scafSeq newBGIseq500_1.fq.gz newBGIseq500_2.fq.gz >newBGIseq500.sam
../../software/samtools view -@ 8 -bS newBGIseq500.sam >newBGIseq500.bam
../../software/samtools sort -@ 8 -O BAM -o  newBGIseq500.sort.bam newBGIseq500.bam
java -jar /home/stu60/software/picard-2.7.1.jar MarkDuplicates I=newBGIseq500.sort.bam O= newBGIseq500.rmdup.bam M= newBGIseq500.rmdup_metrics
../../software/samtools index newBGIseq500.rmdup.bam
/home/stu60/software/bcftools mpileup -O u -f newBGIseq500.scafSeq newBGIseq500.rmdup.bam |/home/stu60/software/bcftools call -Amv -o newBGIseq500.raw.vcf

04.perl
R

data <- read.table("GC_result")
hi <- hist(data$V4, breaks=50, plot=F)
pdf(file="GC.pdf")
plot(hi$mids, hi$density, type="b", col="red", xlab="GC", ylab="Percentage")
dev.off()

Count_scaffold

#!/usr/bin/perl -w
use strict;

open IN, "<newBGIseq500.scafSeq";
my $count=0;
while(<IN>){
    chomp;
    if(/>/){
    $count++;
    }
}
print "The number of scaffolds is : $count\n";
close IN;

GC_content

#!/usr/bin/perl -w
use strict;

open IN,shift;
open OUT,">GC_result";
$/=">";
<IN>;
my $win=100;
my $step=50;
while(<IN>){
    chomp;
    (my $head, my $seq)=split /\n/,$_,2;
    $seq=~s/\n//g;
    my $len=length($seq);
#   print ">$head","\n",$seq,"\n";
    for(my $i=0;$i<=$len-$win;$i+=$step){
        my $GC_rate=0;
        my $win_seq=substr($seq,$i,$win);
        my $GC_Num=($win_seq=~s/G|C/#/g);
        my $N_Num=($win_seq=~s/N/N/g);
        if($N_Num==100){
            $GC_rate=0;
        }
        else{
            $GC_rate=$GC_Num/($win-$N_Num);
        }
        print OUT "$head\t$i\-",$i+$win,"\t$GC_rate\n";
    }

}
close IN;
close OUT;

software

./bzip2 -d bwa-0.7.17.tar.bz2
tar -vf bwa-0.7.17.tar
cd bwa-0.7.17
make