9

The order of -a and -b switched at some point. You want: bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt For reference, this is the end of the help output in version 2.25: Default Output: After each entry in A, reports: 1) The number of features in B that overlapped the A interval. 2) The number of bases in A that ...


6

awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed That produces: chr1 0 10000 LOC101929784::chr1:13201-13800::E1 chr1 10000 10600 LOC101929784::chr1:...


6

If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g. awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \ | mergeBed \ | sed 's/_/\t/' \ | awk -v OFS='\t' '{print $1, $3, $4, $2}' (Assuming "_" is not part of chromosome or identifier name, ...


6

Your setup: import pandas as pd dict1 = {0:['chr1','chr1','chr1','chr1','chr2'], 1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100], 3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'], 4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']} df1 = pd.DataFrame(dict1) print(df1) ## 0 1 2 3 4 5 ## 0 chr1 1 100 ...


5

There are multiple ways to go about it. On the command line you can make a 1 line BED file: chr1 11868 12227 And then bedtools intersect with it. In R, you could load your original BED file and use GenomicRanges: library("GenomicRanges") bed = read.delim("foo.bed", header=F) # Rename this # N.B., BED files use 0-based coordinates, I've switched to 1-...


5

Although you don't mention it, I'm guessing you're using bedtools v2.26.0. Version 2.26.0 of groupBy has a bug in it, which you've encountered (it was fixed shortly after release, so you'll either have to use a version before the bug was introduced, or compile the current source code yourself from https://github.com/arq5x/bedtools2) v2.26.0: local10:~/...


5

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.


4

You could do this with the CGAT toolkit: cgat bed2bed --method=merge --merge-by-name -I bed_with_gene_ids.bed Installing such a massive package might be overkill for this task though.


4

Here is a fun little python script I cooked up for the occasion. It will take a standard fasta file as a command-line argument and turn it into the proper bed format. Using your example fasta: import pandas as pd import sys inFasta = sys.argv[1] # take fasta as command argument def fastaParser(fasta): headers = [] with open(fasta) as f: ...


4

I'm assuming the FASTA header contains all the information you want. You can easily generate a proper bed file with Biopython: from Bio import SeqIO for record in SeqIO.parse('test.fa', 'fasta'): lst = record.id.split(':') chr_, (start, end) = lst[0], lst[1].split('-') print('\t'.join((chr_, start, end))) Note that native Python code is only ...


4

Note: not yet tested, so there may be some additional fiddling with command line options needed The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions): samtools depth -a in1.bam > depth_in1_both.tsv To split this by forward and reverse, you can use an initial pipe through samtools view to exclude or include reverse-...


4

git clone https://github.com/lh3/htsbox cd htsbox && make ./htsbox pileup -cCf ref.fa aln.bam | less -S This output a VCF containing positions covered by at least one read. Tag ADF and ADR give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.


4

I don't think Pandas has this implemented functionality out-of-the-box. Even if it did, solutions not designed specifically for bioinformatics probably rarely handle intervals on different chromosomes correctly unless you split the intervals by chromosome first. Pandas does handle intervals (see docs for the Interval and IntervalIndex classes), but I've ...


4

If you want to find the intersection between ALL your bedfiles, you can try multiIntersectBed (available since bedtools 2.14.3). It should work like this : multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed > answer.bed To sort your ...


4

pybedtools assumes that bedtools is in your path and bedtools itself will return the version with bedtools --version. So: import subprocess subprocess.check_output(['bedtools', '--version'], text=True).strip().split()[1] For me that returns 'v2.30.0'.


3

Here is a fairly simple (and hopefully readable) native Python solution. It assumes that the bed file is sorted prior to parsing: from itertools import groupby class BedRecord(object): def __init__(self, chr, start, end, feature): self.chr = chr self.start = start self.end = end self.feature = feature def parse_records(...


3

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions. Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention: $ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed ...


3

If your headers are always of the form >chrZ:a-b, where chrZ is a chromosome name, a is a zero-indexed, half-open start position, and b is a zero-indexed, half-open stop position: $ awk '($0 ~ /^>/){ print substr($0,2) }' in.fa | sed 's/[:|-]/\t/g' > out.bed If your headers contain coordinates using a different index, then you can make adjustments ...


3

The UCSC format you linked to isn't a BED file, your method should never produce it. What you posted as your desired output is also not a BED file. Below is a BED file: chr1 0 1000 chr1 1000 2000 chr1 2000 3000 chr1 3000 4000 awk and biopython were producing the appropriate thing already. You just need to zoom in on a region in ...


3

You can use bedtools intersect and cut the results: $ cat file1.bed chr1 12048 177033 DUP FALSE -1 22 24 chr1 12048 89237 DUP FALSE 12 10 -1 chr1 2985712 3355300 DEL TRUE 2 2 2 $ cat file2.gtf chr1 22013 41670 gene_id "ENSG00000160075.10" chr1 22013 41670 gene_id "ENSG00000160075.10" chr1 2985732 3355185 ...


3

The following python script will do what you want and should be relatively memory efficient. It processes a single chromosome at a time, so either sort the input or ensure that at least entries in a group are next to each other. #!/usr/bin/env python import sys # Process a sorted list of tuples (e.g., [(0, 10), (10, 20), ...]) # Overlapping items will be ...


3

Via BEDOPS bedops -n and Unix I/O streams: $ echo -e "chr1\t11868\t12227" | bedops -n 1 exon.bed - > answer.bed Or, if you have your genes in a BED file called genes.bed: $ bedops -n 1 exon.bed genes.bed > answer.bed If you have your genes in some other format, like GFF or GTF, you can use gff2bed or gtf2bed, e.g.: $ bedops -n 1 exon.bed <(...


3

Pipe a modified form of the second file to BEDOPS bedmap and the first file, then pipe that result to cut out desired columns: $ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed This should run pretty quickly and use very little memory, as it ...


3

You can do this easily with Hail. Hail primarily uses BED files to annotate genetic datasets (see the last annotate_variants_table example), but you can manipulate BED files using Hail's general facilities for manipulating delimited text files. For example: $ cat genes.bed I 3746 3909 "WBGene00023193" . - I 3746 3909 "WBGene00023193" ...


3

My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues. The linux version of conda installs just fine in WSL and bioconda hae pretty much any bioinformatics package you might need, already ...


3

There is like a thousand different ways how to achieve this. You could use a specialised software for this (like bedtools) or calculate it simply in R. R solution: You can make a function that calculates the number of SNPs in a range and than you can apply it on all the ranges in the table with genomic ranges. snp_table <- read.table('WTSI-OESO_121_1pre....


3

There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.


3

It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the -n option in samtools sort. Otherwise, if you used samtools sort without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered). I think you are ...


3

With v2.27.1 (or any version I ever used) instead of stdin use - (...) | bedtools intersect -a - -b blacklist.bed -v > out.bam Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then ...


Only top voted, non community-wiki answers of a minimum length are eligible