11 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 ...
voiDnyx's user avatar
  • 401
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 ...
winni2k's user avatar
  • 2,256
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 ...
Bioathlete's user avatar
  • 2,574
8 votes

Filtering step for read counts data

If you have your counts in a data.frame called counts, something like this might work: ...
benn's user avatar
  • 3,571
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: ...
terdon's user avatar
  • 9,846
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 ...
Devon Ryan's user avatar
  • 19.6k
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 ...
user172818's user avatar
  • 6,525
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 ...
Devon Ryan's user avatar
  • 19.6k
5 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 ...
terdon's user avatar
  • 9,846
5 votes

PL and QUAL values on VCF file?

The VCF specification provides the definition for the QUAL field. However, QUAL values are often capped by variant callers to a given value. I have seen 100 being used and according to this ...
JRodrigoF's user avatar
  • 827
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: ...
Alex Reynolds's 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]" ...
Pierre's user avatar
  • 1,511
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): ...
Cleb's 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 ...
Devon Ryan's user avatar
  • 19.6k
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.
spvensko's 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 ...
Timur Shtatland's user avatar
4 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: ...
Gawain's user avatar
  • 315
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 ...
Daniel Standage's 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 ...
Kamil S Jaron's user avatar
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 ...
Pierre's user avatar
  • 1,511
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: ...
andremrsantos's 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. ...
Sam Nicholls's 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. ...
story's user avatar
  • 1,573
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 ...
llrs's user avatar
  • 4,693
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? ...
dariober's user avatar
  • 699
3 votes
Accepted

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

You really need to learn some python :) ...
Liam McIntyre's 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: ...
haci's user avatar
  • 4,032
3 votes
Accepted

Filtering file with AWK and writing output to new file

The {} isn't in any way special to the shell or awk. It is used as a replacement by certain other programs, off the top of my ...
terdon's user avatar
  • 9,846
3 votes
Accepted

Filter with bcftools

First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU grep (the default on Linux): ...
terdon's user avatar
  • 9,846
3 votes

dealing with a list of VEP files

For these requirements, it might be possible to simply use filter_vep included in the VEP package for this. Note that multiple --filter flags are ANDed together (i....
Steve's user avatar
  • 3,079

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