Hot answers tagged

18

Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell. If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell. See Csh Programming Considered Harmful and ...


11

TL;DR SH adheres to an official industry standard, but it is not suitable for scientific computing. Bash is considered an informal standard (e.g., by Google). Bash 3 is preferable in most of situations in the world of bioinformatics. Long answer As already described in other answers, SH (/bin/sh, plain Bourne shell, the original UNIX shell) should fully ...


10

Another solution for repeating one command with different parameters is to use parallel. Create a tab-delimited file samples.txt, which contains the sample names and the information for your read groups, e.g.: Sample1 Lib1 Sample2 Lib2 Sample3 Lib3 Then you can use parallel like this: parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:...


8

Something like this might do it: #!/bin/bash # make an array to hold your list of samples/base filenames samples=( S1 S2 S3 S4 S5 ) # make a loop to go through the sample IDs for sample in "${samples[@]}" do : # do something with each sample/file, e.g. make a directory for it mkdir "${sample}"_home # do something else with the sample, e.g. unzip ...


8

I have trialed the free version of InsideDNA, and these were my notes: Cost: \$225/month for a team of 5 with 50TB storage or \$45/month with 10TB storage for individuals (assuming 6 month package: https://insidedna.me/pricing). Software installed: Around 600 bioinformatic tools available and standard command line tools; some popular tools missing (like CD-...


8

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


8

awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa An explanation of the non-obvious bits: BEGIN{FS="::"} split columns using :: as the designator $1~">" If there's a > in the line substr($1,2) Trim the > that was at the beginning of the line.


7

This should work for you: for i in *.tar.gz do dir2="/destination/directory" tar xvzf $i -C $TMPDIR files=($TMPDIR/*) module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam done (If this is the ...


7

The Open Group Base Specifications Issue 7 IEEE Std 1003.1™-2008, 2016 Edition, or "The POSIX Standard" for short, is the standard that defines the interfaces and utilities provided by a Unix system. Among these is the command line shell language and tools (see "Shell & Utilities" in the main index on the page linked above). As far as I know, there is ...


6

First of all, the brace expansion you tried to use isn't valid. While bash can expand {1..3} to 1 2 3, it has no way of knowing how to expand {SRR1002678..SRR1184123}: $ echo {1..3} 1 2 3 $ echo {SRR1002678..SRR1184123} {SRR1002678..SRR1184123} Perhaps you meant SRR{1002678..1184123}? But there's no need for this, you already have the accessions in a file ...


5

if you have a repeating workflow, I strongly recommend to have a look at workflow management systems like snakemake. I've also wrote a little tutorial on biostars about this topic, which might be useful for you. fin swimmer


5

If you just want the number of unique headers, you can do this: grep '>' reference.fasta | sort | uniq | wc -l If you want a list of the unique headers, you can do this: grep '>' reference.fasta | sort | uniq If you want a histogram of how many times each header occurs, you can do this: grep '>' reference.fasta | sort | uniq -c | awk '{printf("%...


5

The generic command sh is quite literally an industry standard, a POSIX standard, to be precise (IEEE 1003.2 and 1003.2a, available for purchase for hundreds of dollars at various websites). In theory, any script that starts with #!/bin/sh should conform to this standard. In practise, most Linux systems have a shell that is close to this standard, but has a ...


5

I would not say bash as a "standard", but it is indeed likely to be the most widely used unix shell and available by default on most modern unix/linux distros. There are a few other more convenient shells like zsh that are broadly compatible with /bin/sh, but they are not as widely available. There is also C-shell and in particular its open-source ...


5

I'm not sure what kinds of bioinformatics tasks you would like to perform, therefore it is difficult to give a good recommendation. If you're specifically working on statistical genetics, I can recommend Hail [1]. Hail is an open-source tool for analyzing genetics data at the tens of terabyte scale. Most of Hail's users do their science in Jupyter notebooks ...


5

Alternatively, you can use the NCBI Entrez Direct UNIX E-utilities Basically, you have to download the install file here: https://www.ncbi.nlm.nih.gov/books/NBK179288/bin/install-edirect.sh In the terminal, install it using: source ./install-edirect.sh Then, you can download your sequence by doing: esearch -db nucleotide -query "NC_030850.1" | efetch -...


5

Assuming all your vcfs are in the same folder: for file in *.vcf; do grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt done


5

Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible. Here's a Perl script that should do exactly ...


5

awk -F "::" '{if($1~">"){gsub(">","");print ">"$2"_"$1} else {print $0}}' foo.fa Basically the same as from Devon but using -F to indicate the initial field delimiter and then using gsub to remove the >.


4

Depending on your applications and uses, you might be interested in checking out CyVerse. It is an NSF funded initiative that provides you with data storage, high performance computing resources, and easy access to commonly used tools. As far as I know, it is free to use once you have an account. I also usually encounter it being used with plant and ...


4

It really depends on what you are trying to do, but here are a few services that I know of. GATK on Google Genomics Cloud: Google and the Broad offer a cloud instance tailored to GATK pipelines. Genomics on Amazon Web Services: I don't think there is anything that makes this unique, but Amazon offers some resources to help get started with genomics/life ...


4

First I installed this libraries: sudo apt-get install libcurl4-openssl-dev libxml2-dev I then installed RCurl: BiocManager::install("RCurl") and lastly: BiocManager::install("DESeq2", version = "3.8")


4

Just wget or curl each as https://www.ncbi.nlm.nih.gov/nuccore/NC_030853.1?report=fasta&log$=seqview&format=text...


4

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"))


4

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f. EDIT AS THE OP HAS EDITED THE QUESTION: First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the ...


4

First of all you are getting an error because you are missing a then in the third line of your code and then a fi to close the if. It should be: for i in *.vcf do if grep -q $i 1.txt; then cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/ fi done And I guess the one-liner below achieves what you are trying to do: while IFS= read -r line; do mv "$line"...


3

The uniq command expects sorted input. Interestingly, the sort command actually has a "unique" option, -u, which means uniq is not strictly needed. For the fastest processing, you can look for the '>' character at the start of lines with grep: grep '^>' reference.fasta | sort -u > reference_headers_unique.fasta For returning the number of unique ...


3

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t"). In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter. egrep -v "/db_xref=\"taxon"| tr "\"" ...


3

Here is a bash script and GNU parallel along with some other tools to split-apply-combine the data. It may or may not work out of the box, but it's not far off and it should give the general idea of how one way of how to approach this. #!/bin/bash # ====== USER PROVIDED PARAMETERS ======== N=20 # number of cores to use ngsfilter="rawdata_scandinavia....


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