8

The GenomicRanges package is a great place to start, something like this: library(GenomicRanges) # make a GenomicRanges set corresponding to the genes genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom") # turn your windows into a data frame and make a GenomicRanges set from that windows = ......


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

In one line, using bedtools zcat Homo_sapiens.GRCh38.93.gtf.gz \ | awk '$3=="gene"' \ | bedtools slop -b 10000 -g contigs.tsv -i - \ | bedtools intersect -u -a intervals.bed -b - This first takes the genes and filters them for the third column being gene. This is important because all GTF lines contain the word gene as gene_id is a manditory attribute ...


5

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits. Convert GTF annotations to a sorted BED file of genes: $ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed Sort your intervals, if unsorted: $ sort-bed intervals.unsorted.bed > intervals....


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

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR. UTR has its own feature in the gtf file. So you can do this: $ awk -v FS="\t" '$3=="UTR"' gencode.gtf If the gtf file is compressed use this instead: $ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"' BTW: Why are you using such an old ...


4

Generically, with BEDOPS: $ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed Or: $ vcf2bed < <(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - > answer.bed Etc. BEDOPS supports input streams; use them where you can. That said, if you're going to query the VCF data frequently, ...


4

this: bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \ java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php java -...


3

First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option -wa, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. ...


3

I was tinkering with the command and was able to complete the execution inverting the order of rsid and bed intervals filtering. The command is as follows: gzcat <myvcf>.vcf.gz | \ java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" | \ java -Xmx10g -jar snpSift.jar intervals <mybed>.bed Maybe my rsid subset is ...


3

This isn't a problem that's easily solved with awk. It's not like you're extracting a feature that's annotated in the GTF file. Instead, you want the empty space between annotated features. A few years ago I wrote a program called LocusPocus for a similar task. It uses a gene annotation to break down a genome into gene loci and intergenic regions. It ...


3

bedtools is my go-to program for operations on genomic intervals. In particular, the bedtools intersect operation is what you're looking for here. $ cat fileA chr1 25 50 chr1 75 200 $ cat fileB chr1 10 60 chr1 80 90 $ bedtools intersect -a fileA -b fileB chr1 25 ...


2

As mentioned by OP, another option is to use pybedtools, which in my opinion is pretty convenient for people already familiar with BedTools. Let's even say df1's format is slightly different than df2: import pandas as pd dict1 = {0: ['chr1', 'chr1', 'chr1', 'chr1', 'chr2'], 1: [1, 100, 150, 900, 1], 2: [100, 200, 500, 950, 100], 3: ['...


2

I implemented pandas "Intervals" and ... it should be a few lines, clearly there are limitations. For non-overlapping data it is very cool however. It will work for overlapping data, BUT if the data you are using as the interval data is overlapping, it falls over. It could work if an independent (non-overlapping) interval was constructed. Anyway, the point ...


2

If you want all transcripts from that gtf file whose type isn't "protein_coding", you can use almost the same command, just change the == ("is") to != ("isn't"): awk '{if($3=="transcript" && $20!="\"protein_coding\";"){print $0}}' gencode.gtf Or, a simpler version: awk '$3=="transcript" && $20!="\"protein_coding\";"' gencode.gtf Note ...


2

This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances. R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're ...


1

The linked gist contains a class I wrote for doing this on GTF files: https://gist.github.com/IanSudery/d8349c22823a475ceb489c3e8aeb448e It uses the GTF class from cgat, which can be found here: https://github.com/cgat-developers/cgat-apps You would use it to do this task like so: from cgat import gtf from cgat import iotools import ...


1

Here is another way to do it in 2019: $ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf Things can be sped up if we compress and index the input.vcf. Then bcftools will have near-immediate random access to the position. $ bgzip -c input.vcf > input.vcf.gz $ tabix input.vcf.gz $ bcftools view -T regions.bed -i 'ID=@<...


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