Gene deserts

Find regions of the mouse genome devoid of any annotation (ESTs, mRNA, repeats, RefSeq and UCSC genes).

Annotation tracks downloaded using the table browser feature of the UCSC Genome Browser. Chromosome sizes of mm9 downloaded from here.

Code for finding regions of 10kb devoid of any annotation.


#!/usr/bin/perl

use strict;
use warnings;

my $bin = 10_000;

my $chrom_file = '/home/davetang/davetang/chrom_size/mm9_chrom_info.txt';
my %chr = ();
open(IN,'<',$chrom_file) || die "Could not open $chrom_file: $!\n";
while(<IN>){
   chomp;
   #chr1    197195432       /gbdb/mm9/mm9.2bit
   my ($chr,$size,$junk) = split(/\s+/);
   $chr{$chr} = $size;
}
close(IN);

my %feature = ();

my @file = qw/mm9_est_gene.bed.gz
mm9_mrna_gene.bed.gz
mm9_refseq_gene.bed.gz
mm9_repeatmasker_gene.bed.gz
mm9_ucsc_gene.bed.gz/;

foreach my $infile (@file){
   open(IN,'-|',"zcat $infile") || die "Could not open $infile: $!\n";
   while(<IN>){
      chomp;
      #chr1    134212712       134221648       BB62599
      my $which_bin = '';
      my ($chr,$start,$end,@junk) = split(/\t/);
      for (my $i = $start; $i <= $end; $i+=$bin){
         $which_bin = int($i / $bin);
         $feature{$chr}{$which_bin} = '1';
      }
      $which_bin = int($end / $bin);
      $feature{$chr}{$which_bin} = '1';
   }
   close(IN);
}

foreach my $chr (keys %chr){
   my $remainder = int($chr{$chr} / $bin);
   foreach (0 .. $remainder){
      next if exists($feature{$chr}{$_});
      my $start = $_ * $bin;
      my $end = $start + 9999;
      my $id = $chr . '_' . $start . '_' . $end;
      print join("\t",$chr,$start,$end,$id),"\n";
   }
}

exit(0);

__END__

In the mm9 genome I found 9,634 10kb regions with no annotation (corresponding to 96,340,000 bp of the genome with no annotation). A few lines of the output (bed format):

chr1 87310000 87319999 chr1_87310000_87319999
chr1 87320000 87329999 chr1_87320000_87329999
chr1 87330000 87339999 chr1_87330000_87339999
chr1 184930000 184939999 chr1_184930000_184939999
chr1 184940000 184949999 chr1_184940000_184949999

Mapping random sized reads to the genome

Simple question, what are the mapping statistics for random pieces of DNA of size 10 to 30 to the human genome?

Generate 1,000,000 random pieces of DNA from size 10 to 30:

#!/usr/bin/perl

use strict;
use warnings;

my @size = (10 .. 30);

foreach my $size (@size){
   warn "$size\n";
   open(OUT,'>',"${size}.fa") || die "Could not open ${size}.fa: $!\n";
   for (1 .. 1000000){
      my $seq = ran_seq($size);
      print OUT ">${_}_$size\n$seq\n";
   }
   close(OUT);
}

sub ran_seq {
   my ($length) = @_;
   my $seq = '';
   my @nuc = qw/A C G T/;
   for (1 .. $length){
      my $rand = int(rand(scalar(@nuc)));
      $seq .= $nuc[$rand];
   }
   return($seq);
}

exit(0);

Map to the hg19 genome using BWA with up to 2 mismatches:

#!/usr/bin/perl

use strict;
use warnings;

my $usage = "Usage: $0  \n";
my $file = shift or die $usage;
my $genome = shift or die $usage;

die "Genome file $genome does not exist\n" unless -e $genome;

my $base = $file;
$base =~ s/\.fa$//;
#my $aln = "bwa aln -n 2 -t 60 $genome $file > $base.sai";
my $aln = "bwa aln -t 60 -n 2 $genome $file > $base.sai";
system($aln);
my $samse = "bwa samse $genome $base.sai $file > $base.sam";
system($samse);
my $bam = "samtools view -S -b $base.sam > $base.bam";
system($bam);
my $bam_sort = "samtools sort $base.bam ${base}_sorted";
system($bam_sort);
my $bam_index = "samtools index ${base}_sorted.bam";
system($bam_index);
unlink("$base.sai");
unlink("$base.sam");
unlink("$base.bam");

system("samstat ${base}_sorted.bam");

exit(0);

__END__

Generate statistics from the BAM file

#!/usr/bin/perl

use strict;
use warnings;

my $usage = "Usage: $0 <infile.bam>\n";
my $infile = shift or die $usage;

my %ed = ();
my %multimappers = ();
my $mapped = '0';
my $unmapped = '0';

open(IN,'-|',"samtools view $infile") || die "Could not open $infile: $!\n";
while(<IN>){
   chomp;
   my ($junk,$flag,@junk) = split(/\t/);
   if ($flag == 4){
      ++$unmapped;
      next;
   } else {
      ++$mapped;
   }
   #965850_13       16      chr1    204577  0       13M     *       0       0       GGCTGCCAGAAGA   *       XT:A:N  NM:i:0  XN:i:13 X0:i:582        XM:i:0  XO:i:0  XG:i:0  MD:Z:13
   if (/NM:i:(\d+)/){
      my $ed = $1;
      if (exists $ed{$ed}){
         $ed{$ed}++;
      } else {
         $ed{$ed} = '1';
      }
   } else {
      die "Could not get mismtach stat on line $.: $_\n";
   }

   if (/X0:i:(\d+)/){
      my $mapped_time = $1;
      if (exists $multimappers{$mapped_time}){
         $multimappers{$mapped_time}++;
      } else {
         $multimappers{$mapped_time} = '1';
      }
   } else {
      die "Could not get multi mapped stat\n";
}
}
close(IN);

my $outfile = $infile;
$outfile =~ s/\.bam/_stat.txt/;
open(OUT,'>',$outfile) || die "Could not write to $outfile: $!\n";

print OUT "Mapped: $mapped\n";
print OUT "Unmapped: $unmapped\n";

print OUT "Edit distances\n";
foreach my $ed (sort {$a <=> $b} keys %ed){
   print OUT "$ed\t$ed{$ed}\n";
}

print OUT "Multimap profile\n";
my %multimap_tally = ();
foreach my $mapped_time (sort {$a <=> $b} keys %multimappers){
   #print OUT "$region\t$region{$region}\n";
   my $count = $multimappers{$mapped_time};
   if ($mapped_time < 11){
      $multimap_tally{$mapped_time} += $count;
   } else {
      $multimap_tally{'42'} += $count;
   }
}

foreach my $mapped (sort {$a <=> $b} keys %multimap_tally){
   my $tally = $multimap_tally{$mapped};
   if ($mapped == 1){
      print OUT "Mapped uniquely\t$tally\n";
   } elsif ($mapped == 42){
      print OUT "Mapped to more than 10 location\t$tally\n";
   } else {
      print OUT "Mapped to $mapped locations\t$tally\n";
   }
}

close(OUT);

exit(0);

Some statistics:

Number of unmapped (2 mismatches) vs. size

10 Unmapped 0
11 Unmapped 1
12 Unmapped 0
13 Unmapped 0
14 Unmapped 0
15 Unmapped 0
16 Unmapped 0
17 Unmapped 6
18 Unmapped 1057
19 Unmapped 16199
20 Unmapped 91842
21 Unmapped 277977
22 Unmapped 545305
23 Unmapped 783820
24 Unmapped 919804
25 Unmapped 974340
26 Unmapped 992416
27 Unmapped 997812
28 Unmapped 999386
29 Unmapped 999828
30 Unmapped 999957

Number of uniquely mapped reads

10 0
11 88
12 5377
13 39204
14 114068
15 163891
16 197747
17 213321
18 230650
19 260335
20 304163
21 345240
22 296730
23 171595
24 70738
25 23848
26 7180
27 2104
28 591
29 166
30 41

Number of perfect matches vs. size

10 1000000
11 999980
12 997445
13 965000
14 847478
15 620991
16 371567
17 162015
18 53240
19 15037
20 3962
21 982
22 255
23 65
24 13
25 5
26 0
27 0
28 0
29 0
30 0

 

The longer the random read, the less likely it can be mapped, which consequently decreases the number of perfect matches. The distribution of uniquely mapped reads peaks at 21 bp reads. Given that after 21 bp the number of unmapped reads increases roughly two-fold, the number of uniquely mapped reads will drop. Thus if a 21 bp sized read is able to map to the genome with 2 mismatches, it is also more likely to be a unique match.

Statistics for individual sizes

cat 10_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 1000000
Multimap profile
Mapped to 5 locations 4
Mapped to 7 locations 4
Mapped to 9 locations 7
Mapped to 10 locations 4
Mapped to more than 11 location 999981

cat 11_sorted_report.txt
Mapped: 999999
Unmapped: 1
Edit distances
0 999980
1 19
Multimap profile
Mapped uniquely 88
Mapped to 2 locations 132
Mapped to 3 locations 254
Mapped to 4 locations 320
Mapped to 5 locations 463
Mapped to 6 locations 516
Mapped to 7 locations 699
Mapped to 8 locations 804
Mapped to 9 locations 973
Mapped to 10 locations 1123
Mapped to more than 11 location 994627

cat 12_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 997445
1 2555
Multimap profile
Mapped uniquely 5377
Mapped to 2 locations 7702
Mapped to 3 locations 9042
Mapped to 4 locations 9714
Mapped to 5 locations 9366
Mapped to 6 locations 8978
Mapped to 7 locations 8264
Mapped to 8 locations 7575
Mapped to 9 locations 6780
Mapped to 10 locations 6388
Mapped to more than 11 location 920814

cat 13_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 965000
1 35000
Multimap profile
Mapped uniquely 39204
Mapped to 2 locations 34571
Mapped to 3 locations 30668
Mapped to 4 locations 29359
Mapped to 5 locations 28932
Mapped to 6 locations 28027
Mapped to 7 locations 27572
Mapped to 8 locations 26393
Mapped to 9 locations 24269
Mapped to 10 locations 22690
Mapped to more than 11 location 708315

cat 14_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 847478
1 152496
2 26
Multimap profile
Mapped uniquely 114068
Mapped to 2 locations 86691
Mapped to 3 locations 65637
Mapped to 4 locations 48070
Mapped to 5 locations 36300
Mapped to 6 locations 27553
Mapped to 7 locations 22235
Mapped to 8 locations 18600
Mapped to 9 locations 16523
Mapped to 10 locations 15280
Mapped to more than 11 location 549043

cat 15_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 620991
1 376676
2 2333
Multimap profile
Mapped uniquely 163891
Mapped to 2 locations 88843
Mapped to 3 locations 61881
Mapped to 4 locations 51557
Mapped to 5 locations 45260
Mapped to 6 locations 41127
Mapped to 7 locations 36060
Mapped to 8 locations 32620
Mapped to 9 locations 28437
Mapped to 10 locations 25541
Mapped to more than 11 location 424783

cat 16_sorted_report.txt
Mapped: 1000000
Unmapped: 0
Edit distances
0 371567
1 600027
2 28406
Multimap profile
Mapped uniquely 197747
Mapped to 2 locations 112470
Mapped to 3 locations 76667
Mapped to 4 locations 54234
Mapped to 5 locations 40978
Mapped to 6 locations 32473
Mapped to 7 locations 26980
Mapped to 8 locations 24282
Mapped to 9 locations 22156
Mapped to 10 locations 20199
Mapped to more than 11 location 391814

cat 17_sorted_report.txt
Mapped: 999994
Unmapped: 6
Edit distances
0 162015
1 704139
2 133834
3 6
Multimap profile
Mapped uniquely 213321
Mapped to 2 locations 110255
Mapped to 3 locations 75728
Mapped to 4 locations 59396
Mapped to 5 locations 48264
Mapped to 6 locations 40014
Mapped to 7 locations 33587
Mapped to 8 locations 28649
Mapped to 9 locations 24931
Mapped to 10 locations 21707
Mapped to more than 11 location 344142

cat 18_sorted_report.txt
Mapped: 998943
Unmapped: 1057
Edit distances
0 53240
1 595190
2 350422
3 91
Multimap profile
Mapped uniquely 230650
Mapped to 2 locations 121417
Mapped to 3 locations 82100
Mapped to 4 locations 62469
Mapped to 5 locations 50773
Mapped to 6 locations 42305
Mapped to 7 locations 36057
Mapped to 8 locations 30848
Mapped to 9 locations 27146
Mapped to 10 locations 23432
Mapped to more than 11 location 291746

cat 19_sorted_report.txt
Mapped: 983801
Unmapped: 16199
Edit distances
0 15037
1 361379
2 606839
3 540
4 6
Multimap profile
Mapped uniquely 260335
Mapped to 2 locations 136022
Mapped to 3 locations 88548
Mapped to 4 locations 63426
Mapped to 5 locations 48358
Mapped to 6 locations 39538
Mapped to 7 locations 33180
Mapped to 8 locations 28925
Mapped to 9 locations 25546
Mapped to 10 locations 22836
Mapped to more than 11 location 237087

cat 20_sorted_report.txt
Mapped: 908158
Unmapped: 91842
Edit distances
0 3962
1 158869
2 743596
3 1705
4 26
Multimap profile
Mapped uniquely 304163
Mapped to 2 locations 149847
Mapped to 3 locations 95531
Mapped to 4 locations 66186
Mapped to 5 locations 49772
Mapped to 6 locations 38053
Mapped to 7 locations 30702
Mapped to 8 locations 24730
Mapped to 9 locations 20523
Mapped to 10 locations 17153
Mapped to more than 11 location 111498

cat 21_sorted_report.txt
Mapped: 722023
Unmapped: 277977
Edit distances
0 982
1 57085
2 661068
3 2847
4 41
Multimap profile
Mapped uniquely 345240
Mapped to 2 locations 145133
Mapped to 3 locations 77801
Mapped to 4 locations 46879
Mapped to 5 locations 29658
Mapped to 6 locations 20092
Mapped to 7 locations 13838
Mapped to 8 locations 9752
Mapped to 9 locations 7108
Mapped to 10 locations 5323
Mapped to more than 11 location 21199

cat 22_sorted_report.txt
Mapped: 454695
Unmapped: 545305
Edit distances
0 255
1 19143
2 432078
3 3165
4 54
Multimap profile
Mapped uniquely 296730
Mapped to 2 locations 87878
Mapped to 3 locations 33860
Mapped to 4 locations 15286
Mapped to 5 locations 7577
Mapped to 6 locations 4182
Mapped to 7 locations 2441
Mapped to 8 locations 1625
Mapped to 9 locations 1057
Mapped to 10 locations 781
Mapped to more than 11 location 3278

cat 23_sorted_report.txt
Mapped: 216180
Unmapped: 783820
Edit distances
0 65
1 6275
2 207910
3 1881
4 49
Multimap profile
Mapped uniquely 171595
Mapped to 2 locations 30212
Mapped to 3 locations 7856
Mapped to 4 locations 2893
Mapped to 5 locations 1318
Mapped to 6 locations 670
Mapped to 7 locations 393
Mapped to 8 locations 255
Mapped to 9 locations 167
Mapped to 10 locations 123
Mapped to more than 11 location 698

cat 24_sorted_report.txt
Mapped: 80196
Unmapped: 919804
Edit distances
0 13
1 1834
2 77582
3 749
4 18
Multimap profile
Mapped uniquely 70738
Mapped to 2 locations 6829
Mapped to 3 locations 1383
Mapped to 4 locations 479
Mapped to 5 locations 211
Mapped to 6 locations 142
Mapped to 7 locations 85
Mapped to 8 locations 66
Mapped to 9 locations 38
Mapped to 10 locations 30
Mapped to more than 11 location 195

cat 25_sorted_report.txt
Mapped: 25660
Unmapped: 974340
Edit distances
0 5
1 562
2 24804
3 283
4 6
Multimap profile
Mapped uniquely 23848
Mapped to 2 locations 1297
Mapped to 3 locations 236
Mapped to 4 locations 97
Mapped to 5 locations 53
Mapped to 6 locations 28
Mapped to 7 locations 22
Mapped to 8 locations 14
Mapped to 9 locations 7
Mapped to 10 locations 12
Mapped to more than 11 location 46

cat 26_sorted_report.txt
Mapped: 7584
Unmapped: 992416
Edit distances
1 136
2 7360
3 86
4 2
Multimap profile
Mapped uniquely 7180
Mapped to 2 locations 283
Mapped to 3 locations 64
Mapped to 4 locations 23
Mapped to 5 locations 13
Mapped to 6 locations 3
Mapped to 7 locations 4
Mapped to 8 locations 1
Mapped to 9 locations 1
Mapped to 10 locations 1
Mapped to more than 11 location 11

cat 27_sorted_report.txt
Mapped: 2188
Unmapped: 997812
Edit distances
1 40
2 2127
3 20
4 1
Multimap profile
Mapped uniquely 2104
Mapped to 2 locations 54
Mapped to 3 locations 13
Mapped to 4 locations 2
Mapped to 5 locations 2
Mapped to 6 locations 2
Mapped to 7 locations 1
Mapped to 8 locations 4
Mapped to 9 locations 1
Mapped to more than 11 location 5

cat 28_sorted_report.txt
Mapped: 614
Unmapped: 999386
Edit distances
1 11
2 598
3 5
Multimap profile
Mapped uniquely 591
Mapped to 2 locations 17
Mapped to 3 locations 2
Mapped to 4 locations 4

cat 29_sorted_report.txt
Mapped: 172
Unmapped: 999828
Edit distances
1 4
2 167
3 1
Multimap profile
Mapped uniquely 166
Mapped to 2 locations 6

cat 30_sorted_report.txt
Mapped: 43
Unmapped: 999957
Edit distances
1 1
2 42
Multimap profile
Mapped uniquely 41
Mapped to 2 locations 2