---
# System prepended metadata

title: 'Project Repo: Analysis on SCARB-2 Gene (XP_021375881)'

---

# Project Repo: Analysis on SCARB-2 Gene (XP_021375881)

Copyright Shane Varghese. Extra Credit Github Repo


## Preamble
SCARB-2 Gene is named "project" for simplicity in editing command lines. All references and version numbers will be located in the Papers Citations. 

*WORK ON ADJUSTING CONTENTS HYPERLINKS*

# 1. Homolog Collection with BLAST

Used NCBI to obtain the BLAST data to collect preliminary data using the accession reference sequence XP_021373098.1

## Commands to Aligning Homologs
Here we downloaded our protein in FASTA format
```
ncbi-acc-download -F fasta -m protein XP_021373098.1
```
Here performed a BLAST search using our protein and are creating a table output for our data
```
blastp -db ../allprotein.fas -query XP_021373098.1.fa -outfmt 0 -max_hsps 1 -out project.blastp.typical.out
blastp -db ../allprotein.fas -query XP_021373098.1.fa -outfmt "6 sseqid pident length mismatch gappen evalue bitscore pident stitle"  -max_hsps 1 -out project.blastp.detail.out
```
Here we used selected our e value cut off. We used an e-value cutoff of 1e-35. We decided to use this e-value to have a high sequence identity, and a high alignment score as well as minimizes false homologs. This results in a filter BLAST sequences
```
awk '{if ($6< 1e-35 )print $1 }' project.blastp.detail.out > project.blastp.detail.filtered.out
```
This was used to count the number of hits after filtering 
```
wc -l project.blastp.detail.filtered.out
```
________________________________________________________________________________________________________________________________________________________________________________
# 2. Sequence Alignment 

Using the filtered BLAST results, we can use MUSCLE to align all the sequences with each other. MUSCLE will create multiple sequence alignments based on different alignment accuracy benchmarks. We used Alv as a visual aid for our alignment. We used alignbuddy to assist with calculating the lengths of the alignments. We used T-Coffee to calculate percent identities among the sequences

## Commands to Align Sequences
We use Seqkit to obtain the BLAST sequences 
```
seqkit grep --pattern-file ~/labs/lab3-$MYGIT/project/project.blastp.detail.filtered.out ~/labs/lab3-$MYGIT/allprotein.fas > ~/labs/lab4-$MYGIT/project/project.homologs.fas
```
This is the muscle command used to create the multiple sequence alignments 
```
muscle -in ~/labs/lab4-$MYGIT/project/project.homologs.fas -out ~/labs/lab4-$MYGIT/project/project.homologs.al.fas
```
Here we use the Alv command to view the alignment and save the alignment in a pdf format 
```
alv -kli --majority ~/labs/lab4-$MYGIT/project/project.homologs.al.fas | less -RS
```
This is how we formatted our alignments on alv. 
```
muscle -in ~/labs/lab4-$MYGIT/project/project.homologs.fas -html -out ~/labs/lab4-$MYGIT/project/project.homologs.al.html
sed -i 's/<PRE>/<pre style="font-family: 'FreeMono', monospaced;">/g' project.homologs.al.html
wkhtmltopdf ~/labs/lab4-$MYGIT/project/project.homologs.al.html --print-media-type ~/labs/lab4-$MYGIT/project/project.homologs.al.pdf
```
This is we begin to use alignbuddy. The First command line provides the width(length) of the alignment. The second command line calculates the alignment after removing any gaps. The last command line calculates the alignment after removing invariant positions.
```
alignbuddy  -al  ~/labs/lab4-$MYGIT/project/project.homologs.al.fas
alignbuddy -trm all  ~/labs/lab4-$MYGIT/gqr/gqr.homologs.al.fas | alignbuddy  -al
alignbuddy -dinv 'ambig' ~/labs/lab4-$MYGIT/gqr/gqr.homologs.al.fas | alignbuddy  -al
```
Here we calculated the average percent identity using the t_coffee program 
```
 alignbuddy -pi ~/labs/lab4-$MYGIT/project/project.homologs.al.fas | awk ' (NR>2)  { for (i=2;i<=NF  ;i++){ sum+=$i;num++} }
     END{ print(100*sum/num) } '
```
________________________________________________________________________________________________________________________________________________________________________________

# 3. Gene Family Sequence Alignment

For the Gene Family Alignment, we used IQ tree to create the optimal phylogenetic trees based on our previously obtained alignments. 

## Commands for the Gene Family Alignments
Here we used IQ-tree with ultra-fast bootstrap support levels using the -bb 1000. This procedure takes about 10 minutes to run.
```
iqtree -s ~/labs/lab5-$MYGIT/project/project.homologs.al.fas -bb 1000 -nt 2 
```
Here we used gotree to create a midpoint rooted tree, since iq tree created an unrooted tree.
```
gotree reroot midpoint -i ~/labs/lab5-$MYGIT/project/project.homologs.al.fas.treefile -o ~/labs/lab5-$MYGIT/project/project.homologs.al.mid.treefile
```
Here we can see the rooted tree using the nw_display program
```
nw_order -c n ~/labs/lab5-$MYGIT/project/project.homologs.al.mid.treefile  | nw_display -
```

# 4. Tree Reconciliation and Rearrangement 

We used the midpoint-rooted gene tree produced in the last section to reconcile the genes using the program Notung. Notung will estimate duplication events and loss events. This allows us to differentiate between orthologs and paralogs. 

## Commands to Reconcile gene and species trees 
This command performs the reconciliation
```
java -jar ~/tools/Notung-3.0-beta/Notung-3.0-beta.jar -s ~/labs/lab5-$MYGIT/species.tre -g ~/labs/lab6-$MYGIT/project/project.homologs.al.mid.treefile --reconcile --speciestag prefix --savepng --events --outputdir ~/labs/lab6-$MYGIT/project/
```
We used this command to assign the internal nodes 
```
grep NOTUNG-SPECIES-TREE ~/labs/lab6-$MYGIT/project/project.homologs.al.mid.treefile.reconciled | sed -e "s/^\[&&NOTUNG-SPECIES-TREE//" -e "s/\]/;/" | nw_display -
```
We used this command to generate a RecPhyloXML to be able to view the gene-within-species using the ThirdKind program
```
python2.7 ~/tools/recPhyloXML/python/NOTUNGtoRecPhyloXML.py -g ~/labs/lab6-$MYGIT/project/project.homologs.al.mid.treefile.reconciled --include.species
thirdkind -Iie -D 40 -f ~/labs/lab6-$MYGIT/project/project.homologs.al.mid.treefile.reconciled.xml -o  ~/labs/lab6-$MYGIT/project/project.homologs.al.mid.treefile.reconciled.svg
```
_________________________________________________________________________________________________________________________________________________________________________

# 5. Protein Domain Prediction

Protein Domain Prediction

## Commands to Aligning Homologs
Here we used RPS-BLAST to identify Pfam Domains without our specific protein sequences. We will be using the project.homolog.fas file containing the original, unaligned sequences from Section 2.
```
sed 's/*//' ~/labs/lab4-$MYGIT/project/project.homologs.fas > ~/labs/lab8-$MYGIT/project/project.homologs.fas
```
Following this, we download the Pfam database from NCBI. Then run the RPS blast command.
```
rpsblast -query ~/labs/lab8-$MYGIT/project/project.homologs.fas -db ~/data/Pfam -out ~/labs/lab8-$MYGIT/project/project.rps-blast.out  -outfmt "6 qseqid qlen qstart qend evalue stitle" -evalue .0000000001
```
Here we used this Rscript to run R without opening an R console. plotTreeandDomain is a program developed by Dr. Rest. We used the data from the tree file and the rps-BLAST output.
```
sudo /usr/local/bin/Rscript  --vanilla ~/labs/lab8-$MYGIT/plotTreeAndDomains.r ~/labs/lab5-$MYGIT/project/project.homologs.al.mid.treefile ~/labs/lab8-$MYGIT/project/project.rps-blast.out ~/labs/lab8-$MYGIT/project/project.tree.rps.pdf
```
We used this command to examine the distribution of Pfam domains across the proteins.
```
cut -f 1 ~/labs/lab8-$MYGIT/project/project.rps-blast.out | sort | uniq -c
cut -f 6 ~/labs/lab8-$MYGIT/project/project.rps-blast.out | sort | uniq -c
```
_______________________________________________________________________________________________________________________________________________________________________________