https://bioinformaticsonline.com/snippets/owner/abhimanyu? https://bioinformaticsonline.com/snippets/view/43735/bash-script-to-extract-sequence-by-ids Tue, 01 Feb 2022 23:20:42 -0600 https://bioinformaticsonline.com/snippets/view/43735/bash-script-to-extract-sequence-by-ids Use a Perl one-liner, grep and seqtk subseq to extract the desired fasta sequences: # Create test input: cat > in.fasta <<EOF >BGI_novel_T016697 Solyc03g033550.3.1 CTGACGTATACAATTAAGCCGCG >BGI_novel_T016313 Solyc03g025570.2.1 TTCAAGTGTTAGTTTCACATCAT >BGI_novel_T018109 Solyc03g080075.1.1 GCAAGGGAAAGAAGTATTACTAG >BGI_novel_T016817 BGI_novel_G001220 GCCCAAGTCATAGGTAGTGCCTG >BGI_novel_T016141 Solyc03g007600.3.1 ACGTACGTACGTACGTACGTACG EOF cat > gene_ids.txt <<EOF Solyc03g033550.3.1 Solyc03g080075.1.1 Solyc00g256710.2.1 Solyc01g010890.3.1 EOF # Extract ids and gene ids into a tsv file: perl -lne '@f = /^>(\S+)\s+(\S+)/ and print join "\t", @f;' in.fasta > ids_gene_ids.tsv # Select ids that correspond to the desired gene ids: grep -f gene_ids.txt ids_gene_ids.tsv | cut -f1 > ids.selected.txt # Extract fasta sequence that correspond to desired gene ids: seqtk subseq in.fasta ids.selected.txt > out.fasta cat out.fasta Output: >BGI_novel_T016697 Solyc03g033550.3.1 CTGACGTATACAATTAAGCCGCG >BGI_novel_T018109 Solyc03g080075.1.1 GCAAGGGAAAGAAGTATTACTAG Note that seqtk can be installed, for example, using conda.]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/43731/install-mgsc-on-ubuntu Mon, 31 Jan 2022 07:50:56 -0600 https://bioinformaticsonline.com/snippets/view/43731/install-mgsc-on-ubuntu (abhiENV) [abhi@hn1 output_OM824]$ conda create -n mgsc python=3.6 graphviz numpy Collecting package metadata (current_repodata.json): done Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source. Collecting package metadata (repodata.json): done Solving environment: done ## Package Plan ## environment location: /home/abhi/anaconda3/envs/mgsc added / updated specs: - graphviz - numpy - python=3.6 The following packages will be downloaded: package | build ---------------------------|----------------- _openmp_mutex-4.5 | 1_gnu 22 KB cairo-1.14.12 | h7636065_2 905 KB fontconfig-2.12.6 | h49f89f6_0 221 KB freetype-2.8 | hab7d2ae_1 542 KB glib-2.69.1 | h4ff587b_1 1.7 MB graphviz-2.40.1 | h25d223c_0 6.5 MB harfbuzz-1.7.6 | h5f0a787_1 474 KB intel-openmp-2022.0.1 | h06a4308_3633 4.2 MB pango-1.42.0 | h377f3fa_0 458 KB pixman-0.40.0 | h7f8727e_1 373 KB readline-8.1.2 | h7f8727e_1 354 KB sqlite-3.37.0 | hc218d9a_0 999 KB wheel-0.37.1 | pyhd3eb1b0_0 33 KB zlib-1.2.11 | h7f8727e_4 108 KB ------------------------------------------------------------ Total: 16.8 MB The following NEW packages will be INSTALLED: _libgcc_mutex pkgs/main/linux-64::_libgcc_mutex-0.1-main _openmp_mutex pkgs/main/linux-64::_openmp_mutex-4.5-1_gnu blas pkgs/main/linux-64::blas-1.0-mkl ca-certificates pkgs/main/linux-64::ca-certificates-2021.10.26-h06a4308_2 cairo pkgs/main/linux-64::cairo-1.14.12-h7636065_2 certifi pkgs/main/linux-64::certifi-2021.5.30-py36h06a4308_0 expat pkgs/main/linux-64::expat-2.4.1-h2531618_2 fontconfig pkgs/main/linux-64::fontconfig-2.12.6-h49f89f6_0 freetype pkgs/main/linux-64::freetype-2.8-hab7d2ae_1 fribidi pkgs/main/linux-64::fribidi-1.0.10-h7b6447c_0 glib pkgs/main/linux-64::glib-2.69.1-h4ff587b_1 graphite2 pkgs/main/linux-64::graphite2-1.3.14-h23475e2_0 graphviz pkgs/main/linux-64::graphviz-2.40.1-h25d223c_0 harfbuzz pkgs/main/linux-64::harfbuzz-1.7.6-h5f0a787_1 icu pkgs/main/linux-64::icu-58.2-he6710b0_3 intel-openmp pkgs/main/linux-64::intel-openmp-2022.0.1-h06a4308_3633 jpeg pkgs/main/linux-64::jpeg-9d-h7f8727e_0 ld_impl_linux-64 pkgs/main/linux-64::ld_impl_linux-64-2.35.1-h7274673_9 libffi pkgs/main/linux-64::libffi-3.3-he6710b0_2 libgcc-ng pkgs/main/linux-64::libgcc-ng-9.3.0-h5101ec6_17 libgomp pkgs/main/linux-64::libgomp-9.3.0-h5101ec6_17 libpng pkgs/main/linux-64::libpng-1.6.37-hbc83047_0 libstdcxx-ng pkgs/main/linux-64::libstdcxx-ng-9.3.0-hd4cf53a_17 libtiff pkgs/main/linux-64::libtiff-4.2.0-h85742a9_0 libtool pkgs/main/linux-64::libtool-2.4.6-h7b6447c_1005 libwebp-base pkgs/main/linux-64::libwebp-base-1.2.0-h27cfd23_0 libxcb pkgs/main/linux-64::libxcb-1.14-h7b6447c_0 libxml2 pkgs/main/linux-64::libxml2-2.9.12-h03d6c58_0 lz4-c pkgs/main/linux-64::lz4-c-1.9.3-h295c915_1 mkl pkgs/main/linux-64::mkl-2020.2-256 mkl-service pkgs/main/linux-64::mkl-service-2.3.0-py36he8ac12f_0 mkl_fft pkgs/main/linux-64::mkl_fft-1.3.0-py36h54f3939_0 mkl_random pkgs/main/linux-64::mkl_random-1.1.1-py36h0573a6f_0 ncurses pkgs/main/linux-64::ncurses-6.3-h7f8727e_2 numpy pkgs/main/linux-64::numpy-1.19.2-py36h54aff64_0 numpy-base pkgs/main/linux-64::numpy-base-1.19.2-py36hfa32c7d_0 openssl pkgs/main/linux-64::openssl-1.1.1m-h7f8727e_0 pango pkgs/main/linux-64::pango-1.42.0-h377f3fa_0 pcre pkgs/main/linux-64::pcre-8.45-h295c915_0 pip pkgs/main/linux-64::pip-21.2.2-py36h06a4308_0 pixman pkgs/main/linux-64::pixman-0.40.0-h7f8727e_1 python pkgs/main/linux-64::python-3.6.13-h12debd9_1 readline pkgs/main/linux-64::readline-8.1.2-h7f8727e_1 setuptools pkgs/main/linux-64::setuptools-58.0.4-py36h06a4308_0 six pkgs/main/noarch::six-1.16.0-pyhd3eb1b0_0 sqlite pkgs/main/linux-64::sqlite-3.37.0-hc218d9a_0 tk pkgs/main/linux-64::tk-8.6.11-h1ccaba5_0 wheel pkgs/main/noarch::wheel-0.37.1-pyhd3eb1b0_0 xz pkgs/main/linux-64::xz-5.2.5-h7b6447c_0 zlib pkgs/main/linux-64::zlib-1.2.11-h7f8727e_4 zstd pkgs/main/linux-64::zstd-1.4.9-haebb681_0 Proceed ([y]/n)? y Downloading and Extracting Packages fontconfig-2.12.6 | 221 KB | ########################################################################################################################################################## | 100% pixman-0.40.0 | 373 KB | ########################################################################################################################################################## | 100% harfbuzz-1.7.6 | 474 KB | ########################################################################################################################################################## | 100% cairo-1.14.12 | 905 KB | ########################################################################################################################################################## | 100% pango-1.42.0 | 458 KB | ########################################################################################################################################################## | 100% sqlite-3.37.0 | 999 KB | ########################################################################################################################################################## | 100% graphviz-2.40.1 | 6.5 MB | ########################################################################################################################################################## | 100% zlib-1.2.11 | 108 KB | ########################################################################################################################################################## | 100% glib-2.69.1 | 1.7 MB | ########################################################################################################################################################## | 100% wheel-0.37.1 | 33 KB | ########################################################################################################################################################## | 100% _openmp_mutex-4.5 | 22 KB | ########################################################################################################################################################## | 100% readline-8.1.2 | 354 KB | ########################################################################################################################################################## | 100% intel-openmp-2022.0. | 4.2 MB | ########################################################################################################################################################## | 100% freetype-2.8 | 542 KB | ########################################################################################################################################################## | 100% Preparing transaction: done Verifying transaction: done Executing transaction: done # # To activate this environment, use # # $ conda activate mgsc # # To deactivate an active environment, use # # $ conda deactivate mgsc -i [path to your assembly graph] -o [output directory name]]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/43730/install-ruby-using-conda Mon, 31 Jan 2022 07:21:48 -0600 https://bioinformaticsonline.com/snippets/view/43730/install-ruby-using-conda [abhi@hn1 MetaCarvel]$ conda install -c conda-forge ruby Collecting package metadata (current_repodata.json): done Solving environment: done ## Package Plan ## environment location: /home/abhi/anaconda3/envs/abhi added / updated specs: - ruby The following packages will be downloaded: package | build ---------------------------|----------------- gdbm-1.18 | h0a1914f_2 190 KB conda-forge ruby-2.6.6 | he592edb_2 5.5 MB conda-forge ------------------------------------------------------------ Total: 5.7 MB The following NEW packages will be INSTALLED: gdbm conda-forge/linux-64::gdbm-1.18-h0a1914f_2 ruby conda-forge/linux-64::ruby-2.6.6-he592edb_2 The following packages will be SUPERSEDED by a higher-priority channel: ca-certificates pkgs/main::ca-certificates-2021.10.26~ --> conda-forge::ca-certificates-2021.10.8-ha878542_0 openssl pkgs/main::openssl-1.1.1m-h7f8727e_0 --> conda-forge::openssl-1.1.1l-h7f98852_0 Proceed ([y]/n)? y Downloading and Extracting Packages gdbm-1.18 | 190 KB | ########################################################################################################################################################## | 100% ruby-2.6.6 | 5.5 MB | ########################################################################################################################################################## | 100% Preparing transaction: done Verifying transaction: done Executing transaction: done]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/39847/convert-fastq-to-fastq Thu, 15 Aug 2019 00:33:31 -0500 https://bioinformaticsonline.com/snippets/view/39847/convert-fastq-to-fastq # Convert FASTQ to FASTA seqtk seq -a IN.fastq > OUT.fasta # Convert FASTQ to FASTA and set bases of quality lower than 20 to N seqtk seq -aQ64 -q20 -n N IN.fastq > OUT.fasta # Download Seqtk https://github.com/lh3/seqtk]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/39687/install-blast-locally Tue, 09 Jul 2019 15:58:10 -0500 https://bioinformaticsonline.com/snippets/view/39687/install-blast-locally Download page https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download Ubuntu / Debian sudo apt-get install ncbi-blast+ CentOS # download latest BLAST version wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.8.1+-x64-linux.tar.gz # decompress tar -zxvf ncbi-blast-2.8.1+-x64-linux.tar.gz # add BLAST location to system PATH export PATH=$HOME/tools/BLAST/ncbi-blast-2.8.1+/bin:$PATH]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/37516/perl-script-to-extract-sequence-by-ids-from-multifasta-file Thu, 09 Aug 2018 10:39:27 -0500 https://bioinformaticsonline.com/snippets/view/37516/perl-script-to-extract-sequence-by-ids-from-multifasta-file #!/usr/bin/perl -w use strict; my $idsfile = "$ARGV[0]"; my $seqfile = "$ARGV[1]"; my %ids = (); open FILE, $idsfile; while(<FILE>) { chomp; $ids{$_} += 1; } close FILE; local $/ = "\n>"; # read by FASTA record open FASTA, $seqfile; while (<FASTA>) { chomp; my $seq = $_; my ($id) = $seq =~ /^>*(\S+)/; # parse ID as first word in FASTA header if (exists($ids{$id})) { $seq =~ s/^>*.+\n//; # remove FASTA header $seq =~ s/\n//g; # remove endlines print ">$id\n$seq\n"; } } close FASTA;]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/36916/perl-script-to-check-fastq-reads-qualities Tue, 12 Jun 2018 04:42:20 -0500 https://bioinformaticsonline.com/snippets/view/36916/perl-script-to-check-fastq-reads-qualities #!/usr/bin/env perl use strict; use warnings; sub readfq { my ($fh, $aux) = @_; @$aux = [undef, 0] if (!defined(@$aux)); return if ($aux->[1]); if (!defined($aux->[0])) { while (<$fh>) { chomp; if (substr($_, 0, 1) eq '>' || substr($_, 0, 1) eq '@') { $aux->[0] = $_; last; } } if (!defined($aux->[0])) { $aux->[1] = 1; return; } } my $name = /^.(\S+)/? $1 : ''; my $seq = ''; my $c; $aux->[0] = undef; while (<$fh>) { chomp; $c = substr($_, 0, 1); last if ($c eq '>' || $c eq '@' || $c eq '+'); $seq .= $_; } $aux->[0] = $_; $aux->[1] = 1 if (!defined($aux->[0])); return ($name, $seq) if ($c ne '+'); my $qual = ''; while (<$fh>) { chomp; $qual .= $_; if (length($qual) >= length($seq)) { $aux->[0] = undef; return ($name, $seq, $qual); } } $aux->[1] = 1; return ($name, $seq); } my @aux = undef; my ($name, $seq, $qual); my ($n, $slen, $qlen) = (0, 0, 0); while (($name, $seq, $qual) = readfq(\*STDIN, \@aux)) { ++$n; $slen += length($seq); $qlen += length($qual) if ($qual); print join("\t", $n, $slen, $qlen), "\n"; } #print join("\t", $n, $slen, $qlen), "\n"; __END__ my @aux = undef; # this is for keeping intermediate data while (my ($name, $seq, $qual) = readfq(\*STDIN, \@aux)) { if( (length($seq) >= 21) && (length($seq) <= 25) ) { print "@$name\n"; print "$seq\n"; print "+\n"; print "$qual\n"; } }]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/36915/perl-script-to-find-palindromic-regions-in-dna-sequences Tue, 12 Jun 2018 04:34:12 -0500 https://bioinformaticsonline.com/snippets/view/36915/perl-script-to-find-palindromic-regions-in-dna-sequences use strict; use warnings; my $pp = qr/(?: (\w) (?1) \g{-1} | \w? )/ix; my $filename = $ARGV[0]; open(my $fh, '<:encoding(UTF-8)', $filename) or die "Could not open file '$filename' $!"; local $/ = ''; while (<$fh>) { chomp; my ($header, @lines) = split "\n"; my $data = join '', @lines; print "$header\n$data\n"; while ($data =~ /(?=($pp))/g) { my $end=($-[0]+length($1)); my $n=(length($1)/2); my $len=length($1); my $midPoint = ($n == int $n) ? $n : int($n + 1); $midPoint=$midPoint+$-[0]; print "$-[0]\t$midPoint\t$end\t$1\t$len\n" if length($1) > 100; } } __DATA__ >TRE|Q47404|Q47404 (409 AA) Glycosyl transferase [Escherichia coli] MIFDASLKKLRKLFVNPIGFFRDSWFFNSKNKAEELLSPLKIKSKNIFIVAHLGQLKKAE LFIQKFSRRSNFLIVLATKKNTEMPRLILEQMNKKLFSSYKLLFIPTEPNTFSLKKVIWF YNVYKYIVLNSKAKDAYFMSYAQHYAIFIWLFKKNNIRCSLIEEGTGTYKTEKKKPLVNI NFYSWIINSIILFHYPDLKFENVYGTFPNLLKEKFDAKKIFEFKTIPLVKSSTRMDNLIH >seq1 TGAATTACTAGAAGTACTTAAAATGATGGTTGGAGGAAATATTCTTGATGATCAAATTGC CGTTAAACTAGGATTTCTTATAAAGGAGGTTGGTAGTAAAATTCATGAAGATCATTAAGT >TRE|Q8VRL9|Q8VRL9 (492 AA) SiaD [Neisseria meningitidis] MLQKIRKALFHPKKFFQDSQWFATPLFSSFAPKSNLFIISTFAQLNQAHSLTKMQKLKNN LLVILYTTQNMKMPKLIQKSVDKELFSVTYMFELPRKPGIVSPKKFLYIQRGYKKLLKTI QPAHLYVMSFAGHYSSLLSLAKKMNITTHLVEEGTATYAPLLESFTYKPTKFEQRFVGNN LHQKGYFDKFDILHVAFPEYAKKIFNANEYHRFFAHSGGISTSQSIAKIQDKYRISQNDY]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/35417/plot-the-density-of-genes-in-r Fri, 02 Feb 2018 03:19:16 -0600 https://bioinformaticsonline.com/snippets/view/35417/plot-the-density-of-genes-in-r #column1 = chromosome name and column2 = start position of the gene # check if ggplot2 is installed, if so, load it, # if not, install and load it if("ggplot2" %in% rownames(installed.packages())){ library(ggplot2) } else { install.packages("ggplot2") library(ggplot2) } # import a text file with gene positions # columns should be: chr, position (no end or gene name required) genes <- read.table("genes.txt",sep="\t",header=T) # make sure the chromosomes are ordered in the way you want # them to appear in the plot genes$chr <- with(genes, factor(chr, levels=paste("chr",c(1:22,"X","Y"),sep=""), ordered=TRUE)) # make a density plot of genes over the provided chromosomes (or scaffolds ...) plottedGenes <- ggplot(genes) + geom_histogram(aes(x=pos),binwidth=1000000) + facet_wrap(~chr,ncol=2) + ggtitle("RefSeq genes density over human genome 19") + xlab("Genomic position (bins 1 Mb)") + ylab("Number of genes") # save it to an image png("genes.png",width=1000,height=1500) print(plottedGenes) dev.off()]]> Abhimanyu Singh https://bioinformaticsonline.com/snippets/view/34233/convert-newline-formated-sequence-into-fasta-format-with-perl Wed, 08 Nov 2017 08:44:06 -0600 https://bioinformaticsonline.com/snippets/view/34233/convert-newline-formated-sequence-into-fasta-format-with-perl use strict; use warnings; my $filename = $ARGV[0]; open(my $fh, '<:encoding(UTF-8)', $filename) or die "Could not open file '$filename' $!"; while (my $row = <$fh>) { chomp $row; my $name= ">seq"."$."; print "$name\n$row\n"; }]]> Abhimanyu Singh