6
==========================================================================
7
7.1 Mask all regions in a genome except for targeted capture regions.
8
==========================================================================
9
# Add 500 bp up and downstream of each probe
11
slopBed -i probes.bed -b 500 > probes.500bp.bed
13
# Get a BED file of all regions not covered by the probes (+500 bp up/down)
15
complementBed -i probes.500bp.bed -g hg18.genome > probes.500bp.complement.bed
17
# Create a masked genome where all bases are masked except for the probes +500bp
19
maskFastaFromBed -in hg18.fa -bed probes.500bp.complement.bed -fo hg18.probecomplement.
23
==========================================================================
24
7.2 Screening for novel SNPs.
25
==========================================================================
26
# Find all SNPs that are not in dbSnp and not in the latest 1000 genomes calls
28
intersectBed -a snp.calls.bed -b dbSnp.bed -v | intersectBed -a stdin -b 1KG.bed
29
-v > snp.calls.novel.bed
33
==========================================================================
34
7.3 Computing the coverage of features that align entirely within an
36
==========================================================================
37
# By default, coverageBed counts any feature in A that overlaps B by >= 1 bp. If
38
you want to require that a feature align entirely within B for it to be counted,
39
you can first use intersectBed with the "-f 1.0" option.
41
intersectBed -a features.bed -b windows.bed -f 1.0 | coverageBed -a stdin -b
42
windows.bed > windows.bed.coverage
45
==========================================================================
46
7.4 Computing the coverage of BAM alignments on exons.
47
==========================================================================
48
# One can combine SAMtools with BEDtools to compute coverage directly from the BAM
49
data by using bamToBed.
51
bamToBed -i reads.bam | coverageBed -a stdin -b exons.bed > exons.bed.coverage
53
# Take it a step further and require that coverage be from properly-paired reads.
55
samtools view -bf 0x2 reads.bam | bamToBed -i stdin | coverageBed -a stdin -b
56
exons.bed > exons.bed.proper.coverage
60
==========================================================================
61
7.5 Computing coverage separately for each strand.
62
==========================================================================
63
# Use grep to only look at forward strand features (i.e. those that end in "+").
65
bamToBed -i reads.bam | grep \+$ | coverageBed -a stdin -b genes.bed >
66
genes.bed.forward.coverage
68
# Use grep to only look at reverse strand features (i.e. those that end in "-").
70
bamToBed -i reads.bam | grep \-$ | coverageBed -a stdin -b genes.bed >
71
genes.bed.forward.coverage
75
==========================================================================
76
7.6 Find structural variant calls that are private to one sample.
77
==========================================================================
80
pairToPair -a sample1.sv.bedpe -b othersamples.sv.bedpe -type neither >
81
sample1.sv.private.bedpe
85
==================================================================================
86
7.7 Exclude SV deletions that appear to be ALU insertions in the reference genome.
87
==================================================================================
88
# We'll require that 90% of the inner span of the deletion be overlapped by a
91
pairToBed -a deletions.sv.bedpe -b ALUs.recent.bed -type notispan -f 0.80 >
92
deletions.notALUsinRef.bedpe
b'\\ No newline at end of file'