|
1、Perl计算N50和序列长度 #/usr/bin/perl -w use strict; 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<@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; } }
#!usr/bin/perl
use warnings;
use strict;
my (@arr, $changdu, %hash, $key);
open (FH, "$ARGV[0]"); #$ARGV[0]指的是脚本运行后跟的第一个参数
while(<FH>){
chomp;
if (/^>/){
$key=$_;
}
$hash{$key}.=$_ unless /^>/;
}
foreach (keys %hash){
$changdu=length($hash{$_});
print "$_\t$changdu\t$hash{$_}\n";
} 2、homo： less -S chr17.vcf | grep -v "^#" | awk '{if($2>=7668402 && $2<=7687550) print $0}' | cut -f10-12 | perl -ne 'chomp;my @a=split/\t/,$_;foreach my $key(@a){my @keys=split/:/,$key;my @keys1=split/\//,$keys[0]; if($keys1[0]==$keys1[1] && $keys1[0]!=0){print "0\t"}elsif($keys1[0]!=$keys1[1]){print "1\t"}else{print "2\t"}}print "\n"' | cut -f1 | awk '{if($1==0) print $0}' | wc –l
Soapdenovo 配置文件如何设置？ 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/stu37/newBGIseq500_1.fq.gz q2=/home/stu37/newBGIseq500_2.fq.gz f1=/home/stu37/newBGIseq500_1.fa.gz f2=/home/stu37/newBGIseq500_2.fa.gz Soapdenovo一步完成： ./
SOAPdenovo all -s lib.cfg -K 35 -D 1 -o ant >>ass.log Fastq转成fasta格式： awk '{if(NR%4 == 1){print
">" substr($0, 2)}}{if(NR%4 == 2){print}}' fastq > fasta 挑取特定突变: grep -v '^#' chr17.vcf |awk '{if($2>=7668402 &&
$2<=7687550){print $0}}' >tp53.vcf vcftools --vcf test.vcf --chr chr17 --from-bp 7668402 --to-bp
7687550 --recode -c > p53.vcf 纯合： less -S tp53.vcf | cut -f 10 | grep
-v "0/0" | awk -F ":" '{print $1}' | awk -F
"\/" '{print $1,$2}' | awk '{if($1==$2) print $0}' | wc –l |
|
|
