Differences

This shows you the differences between two versions of the page.

--- handy_custom_functions [2021/05/12 14:47] – 168.91.18.151
+++ handy_custom_functions [2023/07/25 12:08] (current) – 134.190.232.186
@@ Line 11: / Line 11: @@
 I will discuss here some more custom functions that I found are very useful in my daily workflow. To add these functions to your system, simply add them to your ''.bashrc''
-===Selecting or removing sequences from a FASTA file===
-<code>
+===Reformatting FASTA files downloaded from NCBI===
-# fish out a sequence from a fasta file
-function grabseq {
-    fasta=$3
-    to_grab=$2
-    case "$1" in
-        -s) seqtk subseq $fasta <(grep    "$to_grab" $fasta | sed -e 's/>//' -e 's/ .*//');;
-        -l) seqtk subseq $fasta <(grep -f "$to_grab" $fasta | sed -e 's/>//' -e 's/ .*//');;
-        *) echo -e "grabseq -s <seqname> <fasta> to grab a single sequence\ngrab_seq -l <seqlist> <fasta> to grab a list of sequences"
-    esac
-}
-</code>
-This is essentially a wrapper for the [[https://github.com/lh3/seqtk|seqtk]] tool.
+For most of my analyses, the header format of NCBI FASTA files is very annoying. This function will convert the annoying format into ''>SpeciesName_i_AccessionNumber''
 <code>
-# select any sequence that has <pattern> in the header
+# format NCBI headers to something readable
-# useful if you want to find a single sequence
+function reformat_ncbi_headers {
-$ grabseq -s <pattern> <FASTA>
+    newfile=${1%.fasta}.hdfmt.fasta
+    cp $1 $newfile
-# select a particular set of sequences that have <pattern> in their header and are in the <pattern_list> file
+    sed -i -r -e '/^>/ s/ >.*//' -e 's/>([^ ]*).*\[(.*)\]/>\2_i_\1/' -e 's/ /_/g' $newfile
-# useful if you want to multiple sequences
+    sed -i -r -e '/^>/ s/gi.*ref\|(.*)\|/\1/' $newfile
-$ grabseq -l <pattern_list> <FASTA>
-</code>
-The next function does the opposite of grabseq. It will remove particular sequences from a FASTA file.
-<code>
-# remove a particular entry from a fasta file
-function rmseq {
-    fasta=$3
-    to_rmv=$2
-    case "$1" in
-        -s) seqtk subseq $fasta <( grep ">" $fasta | grep -E -v "$to_rmv" | sed "s/>//" );;
-        -l) seqtk subseq $fasta <( grep ">" $fasta | grep -v -f "$to_rmv" | sed "s/>//" );;
-        *) echo -e "rmseq -s <seqname> <fasta> to remove a single sequence\nremove_seq -l <seqlist> <fasta> to remove a list of sequences"
-    esac
 }
 </code>
-<code>
+===Replacing work names with final names for publication===
-# remove a particular sequence that has <pattern> in the header
-$ rmseq -s <pattern> <FASTA>
-# remove a particular set of sequences that have <pattern> in the header and are in the <pattern_list> file
+In my experience I do my analyses with new genomes / transcriptomes etc I work with 'worknames'. For example 'bin125' or 'L4', or 'BBO'. That is fine while you do your analyses, but in the end when you want to publish your work you'll want to have proper names. This function takes a mapping file that contains the short / worknames in the first column and the final names in the second column, and replaces the worknames with the final names in tree files, FASTA files, you name it.
-$ rmseq -l <pattern_list> <FASTA>
-</code>
-NOTE: Dayana pointed out to me another tool similar to seqtk called [[https://github.com/shenwei356/seqkit|SeqKit]] which has a lot more functionality, essentially making the above custom functions obsolete.
-Selecting sequences:
 <code>
-# select a sequence with the exact <SeqID> as fasta header ID
+# replace taxanames in trees, fasta, etc
-$ seqkit grep -p <SeqID> <FASTA>
+function replace_names {
+    input=$1
-# select a sequence of which the fasta header ID matches a <regex_pattern>
+    mappingfile=$2
-$ seqkit grep -rp <regex_pattern> <FASTA>
+    cp $input $input.nms
+    cat $mappingfile | while read SEARCH REPLACE; do
-# select a set of sequences of which the exact IDs are listed in <SeqID.exact.list>
+        sed -i -r "s/$SEARCH/$REPLACE/" $1.nms
-$ seqkit grep -f <SeqID.exact.list> <FASTA>
+    done
+}
-# select a set of sequences of which the IDs match regex patterns listed in <SeqID.regex.list>
-$ seqkit grep -rf <SeqID.regex.list> <FASTA>
 </code>
-Removing sequences
+===Some other functions===
-<code>
-# remove a sequence with the exact <SeqID> as fasta header ID
-$ seqkit grep -vp <SeqID> <FASTA>
-# remove a sequence of which the fasta header ID matches a <regex_pattern>
-$ seqkit grep -vrp <regex_pattern> <FASTA>
-# remove a set of sequences of which the exact IDs are listed in <SeqID.exact.list>
-$ seqkit grep -vf <SeqID.exact.list>
-# remove a set of sequences of which the IDs match regex patterns listed in <SeqID.regex.list>
-$ seqkit grep -vrf <SeqID.regex.list>
-</code>
-A FASTA header consists of two parts, the header ID and the header DESCRIPTION. The ID is essentially anything between ''>'' and the first space, and the DESCRIPTION is anything after that.
 <code>
-</code>
+# fasta to phylip
+# depends on trimal
+function fa2phy {
+    trimal -in $1 -out ${1%.*}.phylip -phylip
+}
+# reverse complement function
+function revcomp {
+    tr "[ATGCatgcNn]" "[TACGtacgNn]" | rev
+}
-<code>
+# sum up all numbers in a list
-</code>
+function total {
+    tr '\n' '+' $1 | sed "s/\+$/\n/" | bc
+}
-<code>
-</code>
-<code>
 </code>