Hot answers tagged

9

Use perl-rename. This is usually called rename on Debian-based systems like Ubuntu or Mint, and perl-rename or prename on others. Assuming you have it as rename, you can simply do: rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf The -n causes rename to only print what it would do, without doing anything. So once you make sure it does what you ...


6

Perhaps this small test script will help demonstrate some of the principles: #!/usr/bin/env perl use strict; use warnings; use Data::Dumper; my @arr; # define the current line my $line = "foo\tbar\tbaz\n"; # ...


5

just to complete your question, you can do it in R in one line. After setting the directory as your working directory, just type the following: file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf")) Maybe you can try file.rename(list.files(pattern = ".vcf"), ...


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-...


4

Currently you could use either but a major question is which platform will others be using in the future. AFAIK Perl is only superior to Python for regex. Based on the trend I see for new programmers and new software being released: Perl is on the way out and Python is still growing. https://trends.google.com/trends/explore?date=all&q=bioperl,biopython


4

This usually comes down to religious issues, so let me try and steer it back to more objective grounds: What language do you know (better)? Use the library for that one. If you know neither and will be learning a language to use the library, the majority opinion would be that Python is easier to learn. However, some people say that they "click" with Perl ...


4

Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you: Learn the language your colleagues or your advisor use. This way you are able to discuss your code with them and also get help if you run into problems.


4

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop: allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...


4

In the simplest case, if you just want to stop after the first record was printed, you can just add exit (I also corrected the syntax errors you had and added use strict and use warnings; I suggest you get into the habit of using those two, they save you from a lot of grief in the long run): #!/usr/bin/env perl use Bio::DB::GenBank; use Bio::DB::Query::...


4

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP 'TTCG[CGT][ATGC]{1,15}TTCG[CGT]' [CGT] looks for one occurrence of one of the three. [ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases. Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors ...


3

EDIT: (Based on the comments by @terdon, with minor changes) Try this Perl one-liner. It is recommended for general use, even if your file names contain "unexpected" characters such as newlines. perl -e 'for ( 0..$#ARGV ) { rename $ARGV[$_] => "A@{[ $_ + 1 ]}.vcf"; } ' *pass.vcf Here, -e tells perl to look for code on the command line instead of a ...


3

If your filenames do not contain blank spaces, you can do it with a for-loop in bash: Make a test run with echo so that the command is only printed but not yet executed: for i in *.vcf; do echo mv $i A${k}.vcf ; let k++ ; done And if the result looks good remove the echo statement to execute the mv commands: for i in *.vcf; do mv $i A${k}.vcf ; let k++ ; ...


3

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk: $ awk 'NR==FNR{a[$1]=0; next} { if($2 in a){ a[$2]=1 } }END{ for(i in a){ print i,a[i] } }' pangenome.txt blast....


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

You don't need a loop. You can do the whole thing with a simple awk one-liner: awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta I ran this on your files and got: $ for file in gene*; do echo "=== File: $file ==="; cat $file; done === File: gene1.fasta === >locus1 ATGCGTAGAG >ysr ATGTAGCGA >siv TAGTAGTAT === ...


2

A possible approach in plain python: #!/usr/bin/env python3 import sys def both_zero(freq1, freq2): return freq1 == "0" and freq2 == "0" with open(sys.argv[1], "r") as snp_file: for line1 in snp_file: # Consume one extra line to get the frequencies # for the alternative allele line2 = snp_file.readline() # "zip" ...


2

The following Perl documentation pages should be informative: split - for splitting a scalar at all matches of a defined pattern perlreftut - discusses the approach of anonymous variables and how to combine them


2

Here is a way to do it using a python dictionary: #!/usr/bin/env python3 import sys gene_filename = sys.argv[1] blast_filename = sys.argv[2] # Set to 0 all genes in the gene file genes = {} with open(gene_filename, "r") as gene_file: for line in gene_file: genes[line.strip()] = 0 # Set to 1 those genes found in the blast file with open(...


2

Using python it should be relatively trivial, using sets. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the union method of python's sets. pangenome = {x.rstrip() for x in open('pangenome.txt')} blast = {x.rstrip().split()[1] for x in open('blast.txt')} for x in pangenome: ...


2

Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl: #!/usr/bin/perl use strict; use warnings; ## This will hold the information for each strain my %strains; ## These are the functions that will test each strain sub isMDR{...


2

While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks. The quick fix - in python looks like this: #This is the main function which uses Ensembl REST API to access the Enesembl Compara database #And retrive the alignment blocks which correspond to the MM10 standard....


2

Not exactly what you were looking for but may be you can use something like this? ncbi2=NC_002021 efetch -db protein -format gpc -id "$ncbi2" \ | xtract -pattern INSDSeq_feature-table \ -group INSDFeature \ -if INSDFeature_key -equals 'mat_peptide' \ -tab '\n' \ -element INSDFeature_location,INSDQualifier_value \ ...


2

I fear this entire thread is going to veer dangerously close to opinionation. But regardless of how one feels about the Perl language and its various libraries and communities, I think we can objectively say that the era of Perl's dominance in the bioinformatics community has passed. A lot of tools and libraries implemented in Perl are still widely used, and ...


2

I believe a pandas approach would be fast enough (probably faster than any other Python approach without pandas). After converting your table to a pandas data frame, you can do something like df.filter(regex=(".*Chr01")) for the first chromosome then write the resulting data to file with the to_csv() method. import pandas as pd data = pd.read_csv('input....


2

If you can use awk, you could use this simple script: awk ' BEGIN {p=0} NR==FNR {gid[$1];next} /^>/{ p=0; for(g in gid){ re="\\sgene:"g"\\s" if($0~re){p=1;break} } }p' gene_list cds.fasta You can also put the code inside a script file: extract_genes.awk #!/usr/bin/awk BEGIN { p = 0; # p indicates whether to print a line or not } # ...


2

I would align your sequences with the software available here: https://github.com/veg/hyphy-analyses/tree/master/codon-msa


2

If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below: $ cat accs.txt ATO98108.1 ATO98120.1 ATO98132.1 AVP78031.1 AVP78042.1 $ cat accs.txt | epost -db nuccore | efetch -format acc ## no output because none of them are nucleotide accessions $ cat accs.txt | epost -db protein -format acc | efetch -format acc ...


2

Similar to reading files in other languages, it is necessary to trim the newlines from each line. If you add chomp; to the first line of the while loop, the code runs just fine (note that I also print a newline at the end of the loop). open (IDs, '<', "genes.txt") or die "can't read input file"; while (<IDs>){ chomp; my $...


1

Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV. Update, I can provide an answer. It's HDF5 with pandas dataframes (Python). Simply df.to_hdf('path/file.h5', 'key', format='table') This isn't SQL but it works: What it gives, it's fast, that's it's purpose (not my main reason though) it ...


1

If you have the sequences in a seperate file like this: $ cat input LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF you can use a simple cut to ...


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