Skip to content
Snippets Groups Projects
Commit 437a9c80 authored by Aaron's avatar Aaron
Browse files

choose null for mapBed

parent 7826cc4c
No related branches found
No related tags found
No related merge requests found
Showing
with 368 additions and 5 deletions
......@@ -328,7 +328,7 @@ void BedIntersect::IntersectBam(string bamFile) {
BED bed;
MakeBedFromBam(bam, chrom, bed_blocks, bed);
bool overlapsFound = false;
if ((_bamOutput == true) && (_obeySplits == false))
if ((_bamOutput == true) && (_obeySplits == false))
{
overlapsFound = _bedB->anyHits(bed.chrom, bed.start, bed.end,
bed.strand, _sameStrand, _diffStrand,
......
......@@ -19,7 +19,8 @@ double GetUserColumn(const string s);
BedMap::BedMap(string bedAFile, string bedBFile, int column, string operation,
float overlapFraction, bool sameStrand,
bool diffStrand, bool reciprocal,
string nullValue, bool printHeader)
bool choseNullValue, string nullValue,
bool printHeader)
{
_bedAFile = bedAFile;
......@@ -33,6 +34,8 @@ BedMap::BedMap(string bedAFile, string bedBFile, int column, string operation,
_nullValue = nullValue;
_printHeader = printHeader;
if (!choseNullValue && operation == "count")
_nullValue = "0";
Map();
}
......@@ -84,7 +87,7 @@ string BedMap::MapHits(const BED &a, const vector<BED> &hits) {
output << vo.GetMode();
else if (_operation == "antimode")
output << vo.GetAntiMode();
else if (_operation == "count")
else if (_operation == "count")
output << setprecision (PRECISION) << vo.GetCount();
else if (_operation == "count_distinct")
output << setprecision (PRECISION) << vo.GetCountDistinct();
......
......@@ -41,7 +41,8 @@ public:
BedMap(string bedAFile, string bedBFile, int column, string operation,
float overlapFraction, bool sameStrand,
bool diffStrand, bool reciprocal,
string nullValue, bool printHeader);
bool choseNullValue, string nullValue,
bool printHeader);
// destructor
~BedMap(void);
......
......@@ -48,6 +48,7 @@ int map_main(int argc, char* argv[]) {
bool sameStrand = false;
bool diffStrand = false;
bool printHeader = false;
bool choseNullValue = false;
// check to see if we should print out some help
if(argc <= 1) showHelp = true;
......@@ -114,6 +115,7 @@ int map_main(int argc, char* argv[]) {
}
else if (PARAMETER_CHECK("-null", 5, parameterLength)) {
nullValue = argv[i + 1];
choseNullValue = true;
i++;
}
else if(PARAMETER_CHECK("-header", 7, parameterLength)) {
......@@ -146,7 +148,8 @@ int map_main(int argc, char* argv[]) {
BedMap *bm = new BedMap(bedAFile, bedBFile, column, operation,
overlapFraction, sameStrand,
diffStrand, reciprocalFraction,
nullValue, printHeader);
choseNullValue, nullValue,
printHeader);
delete bm;
return 0;
}
......
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:249250621
one_blocks 16 chr1 1 40 30M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50
BT=../../bin/bedtools
check()
{
if diff $1 $2; then
echo ok
else
echo fail
fi
}
###########################################################
###########################################################
# BAM files #
###########################################################
###########################################################
samtools view -Sb one_block.sam > one_block.bam 2>/dev/null
samtools view -Sb two_blocks.sam > two_blocks.bam 2>/dev/null
samtools view -Sb three_blocks.sam > three_blocks.bam 2>/dev/null
samtools view -Sb sam-w-del.sam > sam-w-del.bam 2>/dev/null
##################################################################
# Test one block without -split
##################################################################
echo " bamtobed.t1...\c"
echo \
"chr1 0 30 one_blocks 40 -" > exp
$BT bamtobed -i one_block.bam > obs
check obs exp
rm obs exp
##################################################################
# Test one block with -split
##################################################################
echo " bamtobed.t2...\c"
echo \
"chr1 0 30 one_blocks 40 -" > exp
$BT bamtobed -i one_block.bam -split > obs
check obs exp
rm obs exp
##################################################################
# Test two blocks without -split
##################################################################
echo " bamtobed.t3...\c"
echo \
"chr1 0 40 two_blocks 40 -" > exp
$BT bamtobed -i two_blocks.bam > obs
check obs exp
rm obs exp
##################################################################
# Test two blocks with -split
##################################################################
echo " bamtobed.t4...\c"
echo \
"chr1 0 15 two_blocks 40 -
chr1 25 40 two_blocks 40 -" > exp
$BT bamtobed -i two_blocks.bam -split > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks without -split
##################################################################
echo " bamtobed.t5...\c"
echo \
"chr1 0 50 three_blocks 40 -" > exp
$BT bamtobed -i three_blocks.bam > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks with -split
##################################################################
echo " bamtobed.t6...\c"
echo \
"chr1 0 10 three_blocks 40 -
chr1 20 30 three_blocks 40 -
chr1 40 50 three_blocks 40 -" > exp
$BT bamtobed -i three_blocks.bam -split > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks with -bed12
##################################################################
echo " bamtobed.t7...\c"
echo \
"chr1 0 50 three_blocks 40 - 0 50 255,0,0 3 10,10,10 0,20,40" > exp
$BT bamtobed -i three_blocks.bam -bed12 > obs
check obs exp
rm obs exp
##################################################################
# Ensure that both ways of getting blocks from a spliced alignment
# are indenticsl
##################################################################
echo " bamtobed.t8...\c"
$BT bamtobed -i three_blocks.bam -split > split
$BT bamtobed -i three_blocks.bam -bed12 | $BT bed12tobed6 > bed12
check split bed12
rm split bed12
rm *.bam
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:249250621
three_blocks 16 chr1 1 40 10M10N10M10N10M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:249250621
two_blocks 16 chr1 1 40 15M10N15M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50
\ No newline at end of file
chr1 0 50 bed6 1 +
\ No newline at end of file
chr1 0 50 one_blocks_match 0 + 0 0 0 1 50, 0,
BT=../../bin/bedtools
check()
{
if diff $1 $2; then
echo ok
else
echo fail
fi
}
##################################################################
# Test one block
##################################################################
echo " bed12tobed6.t1...\c"
echo \
"chr1 0 50 one_blocks_match 0 +" > exp
$BT bed12tobed6 -i one_blocks.bed > obs
check obs exp
rm obs exp
##################################################################
# Test two blocks
##################################################################
echo " bed12tobed6.t2...\c"
echo \
"chr1 0 10 two_blocks_match 0 +
chr1 40 50 two_blocks_match 0 +" > exp
$BT bed12tobed6 -i two_blocks.bed > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks
##################################################################
echo " bed12tobed6.t3...\c"
echo \
"chr1 0 10 three_blocks_match 0 +
chr1 20 30 three_blocks_match 0 +
chr1 40 50 three_blocks_match 0 +" > exp
$BT bed12tobed6 -i three_blocks.bed > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks and add block numbers
##################################################################
echo " bed12tobed6.t4...\c"
echo \
"chr1 0 10 three_blocks_match 1 +
chr1 20 30 three_blocks_match 2 +
chr1 40 50 three_blocks_match 3 +" > exp
$BT bed12tobed6 -i three_blocks.bed -n > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks and add block numbers. Test reverse strand
##################################################################
echo " bed12tobed6.t5...\c"
echo \
"chr1 0 10 three_blocks_match 3 -
chr1 20 30 three_blocks_match 2 -
chr1 40 50 three_blocks_match 1 -" > exp
sed -e 's/\+/\-/' three_blocks.bed | $BT bed12tobed6 -n > obs
check obs exp
rm obs exp
\ No newline at end of file
chr1 0 50 three_blocks_match 0 + 0 0 0 3 10,10,10, 0,20,40,
chr1 0 50 two_blocks_match 0 + 0 0 0 2 10,10, 0,40,
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:1000
one_blocks 16 chr1 1 40 30M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50
\ No newline at end of file
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:1000
one_bp_del 16 chr1 1 40 10M1D10M * 0 0 GAAGGCCACCGCCGCGCCGC CCCDDB?=FJIIJIGIIJIG MD:Z:50
BT=../../bin/bedtools
check()
{
if diff $1 $2; then
echo ok
else
echo fail
fi
}
###########################################################
###########################################################
# BAM files #
###########################################################
###########################################################
samtools view -Sb one_block.sam > one_block.bam 2>/dev/null
samtools view -Sb two_blocks.sam > two_blocks.bam 2>/dev/null
samtools view -Sb three_blocks.sam > three_blocks.bam 2>/dev/null
samtools view -Sb sam-w-del.sam > sam-w-del.bam 2>/dev/null
##################################################################
# Test three blocks without -split
##################################################################
echo " genomecov.t1...\c"
echo \
"chr1 0 50 1" > exp
$BT genomecov -ibam three_blocks.bam -bg > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks with -split
##################################################################
echo " genomecov.t2...\c"
echo \
"chr1 0 10 1
chr1 20 30 1
chr1 40 50 1" > exp
$BT genomecov -ibam three_blocks.bam -bg -split > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks with -split and -bga
##################################################################
echo " genomecov.t3...\c"
echo \
"chr1 0 10 1
chr1 10 20 0
chr1 20 30 1
chr1 30 40 0
chr1 40 50 1
chr1 50 1000 0" > exp
$BT genomecov -ibam three_blocks.bam -bga -split > obs
check obs exp
rm obs exp
##################################################################
# Test blocked BAM from multiple files w/ -bga and w/o -split
##################################################################
echo " genomecov.t4...\c"
echo \
"chr1 0 30 3
chr1 30 40 2
chr1 40 50 1
chr1 50 1000 0" > exp
samtools merge -f /dev/stdout *block*.bam | $BT genomecov -ibam - -bga > obs
check obs exp
rm obs exp
##################################################################
# Test blocked BAM from multiple files w/ -bga and w -split
##################################################################
echo " genomecov.t5...\c"
echo \
"chr1 0 10 3
chr1 10 15 2
chr1 15 20 1
chr1 20 25 2
chr1 25 30 3
chr1 30 50 1
chr1 50 1000 0" > exp
samtools merge -f /dev/stdout *block*.bam | $BT genomecov -ibam - -bga -split > obs
check obs exp
rm obs exp
##################################################################
# Test three blocks with -split and -dz
##################################################################
echo " genomecov.t6...\c"
echo \
"chr1 0 1
chr1 1 1
chr1 2 1
chr1 3 1
chr1 4 1
chr1 5 1
chr1 6 1
chr1 7 1
chr1 8 1
chr1 9 1
chr1 20 1
chr1 21 1
chr1 22 1
chr1 23 1
chr1 24 1
chr1 25 1
chr1 26 1
chr1 27 1
chr1 28 1
chr1 29 1
chr1 40 1
chr1 41 1
chr1 42 1
chr1 43 1
chr1 44 1
chr1 45 1
chr1 46 1
chr1 47 1
chr1 48 1
chr1 49 1" > exp
$BT genomecov -ibam three_blocks.bam -dz -split > obs
check obs exp
rm obs exp
##################################################################
# Test SAM with 1bp D operator
##################################################################
echo " genomecov.t7...\c"
echo \
"chr1 0 10 1
chr1 11 21 1" > exp
$BT genomecov -ibam sam-w-del.bam -bg > obs
check obs exp
rm obs exp
rm *.bam
\ No newline at end of file
@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:1000
three_blocks 16 chr1 1 40 10M10N10M10N10M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50
chr1 0 50 three_blocks_match 0 + 0 0 0 3 10,10,10, 0,20,40,
chr1 10 60 three_blocks_nomatch 0 + 0 0 0 3 11,10,10, 0,20,40,
chr1 10 60 three_blocks_nomatch 0 + 0 0 0 3 10,10,10, 0,20,40,
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment