10 votes
Accepted

Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

I'm trying to figure out why. This will be a longer read, Tl;dr at the end: Doing a match There is a good match of Q against ...
user avatar
  • 391
9 votes

Delete all 4 lines of a fastq read from a fastq file using read ID

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ. For a detailed description of what can go wrong in FASTQ parsing see for ...
user avatar
  • 2,096
8 votes
Accepted

How do I find identical sequences in a FASTA file?

The trick would be to swap the key in the dictionary to be the sequence itself. Also I would recommend using a different separator that "_" since that is what the current ids have so that you can ...
user avatar
  • 2,576
8 votes

Filtering step for read counts data

If you have your counts in a data.frame called counts, something like this might work: ...
user avatar
  • 3,551
6 votes

Delete all 4 lines of a fastq read from a fastq file using read ID

If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command: ...
user avatar
  • 8,161
6 votes
Accepted

How can longest isoforms (per gene) be extracted from a FASTA file?

While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory ...
user avatar
  • 19.2k
6 votes

Fast processing of fastq data

Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose. C++, Rust, Go, Nim and Julia can ...
user avatar
  • 5,764
5 votes
Accepted

How to safely and efficiently convert subset of bam to fastq?

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you ...
user avatar
  • 19.2k
4 votes

Any fast options to query large VCF bed intervals?

Generically, with BEDOPS: $ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed Or: ...
user avatar
4 votes

Any fast options to query large VCF bed intervals?

this: bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \ java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" ...
user avatar
  • 1,473
4 votes

How do I find identical sequences in a FASTA file?

Just as a small variation to @Bioathlete's answer in case you want to write the fasta using Biopython (e.g. to add names and description): ...
user avatar
  • 743
4 votes
Accepted

Select top 100 genes ranked by variance in read counts

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it ...
user avatar
  • 19.2k
4 votes
Accepted

Output from vcftools missingness

Humans are diploid, so you can expect to see up to 2*N (2*68=136) alleles, so N_DATA is the number of observances of that allele.
user avatar
4 votes
Accepted

Edit all the fasta headers using awk

EDIT after OP's update: Try this Perl one-liner: perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file Here, the ...
user avatar
4 votes
Accepted

What are values of FILTER column of vcf files produced by Mutect2

You haven't told us what commands you ran, so I am assuming you first ran Mutect2 and then FilterMutectCalls. If so, the ...
user avatar
  • 8,161
3 votes

Delete all 4 lines of a fastq read from a fastq file using read ID

When I can I like to do file processing line-by-line in the spirit of UNIX tools. You can read 4 lines from a Fastq file into 4 tab-separated values on a single line using ...
user avatar
3 votes

Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern

For at least one in a pair matching the sequence, if the read names of pairs are identical, some grep magic could do the job: Print the sequence at the beginning of the line and read the name in the ...
user avatar
  • 5,327
3 votes

Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern

First option My solution uses the following steps: use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C ) use jjs (java ...
user avatar
  • 1,473
3 votes

Any fast options to query large VCF bed intervals?

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: ...
user avatar
3 votes

How can longest isoforms (per gene) be extracted from a FASTA file?

Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your FASTA so use it. ...
user avatar
3 votes

How can longest isoforms (per gene) be extracted from a FASTA file?

Here is a solution in R. Could get really slow with big files. Works for the example you posted. ...
user avatar
  • 1,498
3 votes

Filtering step for read counts data

While this answers explains how to do it I want to address when and why and which thresholds to do it. Filtering the genes with low counts is usually done because the counts are not reliable it would ...
user avatar
  • 4,622
3 votes

Edit all the fasta headers using awk

Your own answer doesn't seem to do what you ask as it doesn't add the sequential identifier to each sequence. Maybe try this? ...
user avatar
  • 649
3 votes
Accepted

How to filter out partially similar strings from two lists and make one list?

You really need to learn some python :) ...
user avatar
3 votes

How to copy only certain counts for genes in tsv file to new file in linux

Something like this should do job: ...
user avatar
  • 3,477
3 votes
Accepted

How to subset a GRanges object based on a specific genomic window of interest?

Specifying the range of interest (+/- 1 Mb of start coordinate) as a new GRanges object and then subsetting the original GRanges object based on the overlap of genes with the query works: ...
user avatar
  • 115
2 votes

How to safely and efficiently convert subset of bam to fastq?

Ok, I wrote a bit brute pysam / BioPython parser that uses index of bam to get proper order of read pair for R1 / R2 files and bitwise flag. It should not be too difficult to add more sophisticated ...
user avatar
  • 5,327
2 votes

How do I find identical sequences in a FASTA file?

using datamash ans some awk ...
user avatar
  • 1,473
2 votes

Filter out outliers of the scRNA-seq (heterogenous cells)

From what I known, there is no clear consensus in the field and it depends on the type of cells you are interrogating. However, if mitoRNA/endogenousRNA ratio does not fit for your purposes, other ...
user avatar
  • 1,022
2 votes

Edit all the fasta headers using awk

sed -e 's:^\(.*\)_.*_\(.*\)$:\1/\2:' < input > output
user avatar
  • 364

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