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 ...
- 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 ...
- 2,236
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 ...
- 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:
...
- 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:
...
- 9,662
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 ...
- 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 ...
- 6,485
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 ...
- 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 ...
- 9,662
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 ...
- 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:
...
- 3,125
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]"
...
- 1,501
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):
...
- 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 ...
- 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.
- 56
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 ...
- 839
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 ...
- 5,060
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 ...
- 5,517
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 ...
- 1,501
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:
...
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.
...
- 772
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.
...
- 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 ...
- 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?
...
- 689
3
votes
Accepted
How to filter out partially similar strings from two lists and make one list?
You really need to learn some python :)
...
- 579
3
votes
How to copy only certain counts for genes in tsv file to new file in linux
Something like this should do job:
...
- 3,947
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:
...
- 115
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 ...
- 9,662
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):
...
- 9,662
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....
- 3,059
Only top scored, non community-wiki answers of a minimum length are eligible