---
# System prepended metadata

title: PRC1 pipeline
tags: [bio312]

---

# PRC1 pipeline
PCGF1
https://www.ncbi.nlm.nih.gov/gene/5517480

## lab3
```
#cd ~/labs
mkdir ~/prc1
cd ~/prc1
git clone https://github.com/Bio312/lab3
cd lab3
gunzip proteomes/*.gz
cat  proteomes/* > allprotein.fas
makeblastdb -in allprotein.fas -dbtype prot
ncbi-acc-download -F fasta -m protein XP_032217712.2
blastp -db allprotein.fas -query XP_032217712.2.fa -outfmt 0 -max_hsps 1 -out prc1.blastp.typical.out
blastp -db allprotein.fas -query XP_032217712.2.fa -outfmt "6 sseqid pident length mismatch gapopen evalue bitscore pident stitle"  -max_hsps 1 -out prc1.blastp.detail.out
grep Nvectensis prc1.blastp.detail.out | less -S
awk '{if ($6<0.000002)print $1 }' prc1.blastp.detail.out > prc1.blastp.detail.filtered.out
wc -l prc1.blastp.detail.filtered.out
 seqkit grep --pattern-file prc1.blastp.detail.filtered.out allprotein.fas > prc1.homologs.fas
muscle -in prc1.homologs.fas -out prc1.homologs.al.fas
 t_coffee -other_pg seq_reformat -in prc1.homologs.al.fas -output sim  |  less -S
total 12.7% identity!!!

~/tools/mafft-linux64/mafft.bat --localpair --maxiterate 100  prc1.homologs.fas   > prc1.homologs.mafftal.fas 
 t_coffee -other_pg seq_reformat -in prc1.homologs.mafftal.fas -output sim  |  less -S
total 20.64% identity!

alv -kli --majority prc1.homologs.mafftal.fas | less -RS
```
## canonical visual opsins
XP_021357501.1
```
mkdir ~/opsin
cd ~/opsin
    git clone https://github.com/Bio312/lab3
cd lab3
gunzip proteomes/*.gz
cat  proteomes/* > allprotein.fas
makeblastdb -in allprotein.fas -dbtype prot

ncbi-acc-download -F fasta -m protein XP_021357501.1
blastp -db allprotein.fas -query XP_021357501.1.fa -outfmt 0 -max_hsps 1 -out opsin.blastp.typical.out


blastp -db allprotein.fas -query XP_021357501.1.fa -outfmt "6 sseqid pident length mismatch gapopen evalue bitscore pident stitle"  -max_hsps 1 -out opsin.blastp.detail.out


awk '{if ($6< 6e-37 )print $1 }' opsin.blastp.detail.out > opsin.blastp.detail.filtered.out

wc -l opsin.blastp.detail.filtered.out #59

grep -o -E "^[A-Z][a-z]+\." opsin.blastp.detail.filtered.out  | sort | uniq -c

   5 Aplanci.
   6 Avaga.
  14 Bbelcheri.
   2 Cintestinalis.
   7 Dmelanogaster.
   4 Hsapiens.
  10 Lanatina.
  11 Myessoensis.

#expecting from below search...
                                               qseqid tophits    avPid avLength
2:     Myessoensis.XP_021357501.1_rhodopsin_GQcoupled      59 35.38475 307.5424
   Myessoensis Lanatina Bbelcheri Aplanci Cintestinalis Avaga Dmelanogaster
2:          11       10        14       5             1     6             7
   Hsapiens Egranulosus minSp maxSp numZero unchar
2:        5           0     0    14       1  FALSE


grep -o -E "^[A-Z][a-z]+\." opsin.blastp.detail.filtered.out  | sort | uniq -c | wc -l 

8



 seqkit grep --pattern-file opsin.blastp.detail.filtered.out allprotein.fas > opsin.homologs.fas
muscle -in opsin.homologs.fas -out opsin.homologs.al.fas
 t_coffee -other_pg seq_reformat -in opsin.homologs.al.fas -output sim 

25.10% ID

alv -kli --majority opsin.homologs.al.fas | less -RS

```
## lab4
```
iqtree -s opsin.homologs.al.fas  -m VT+F+R10
```
start 2:14, end 2:28; 14 minutes
44 dups, 54 losses (59 genes) score 120 (2.03 points per gene)

## XP_021373098.1 as query (discovered from below)

ncbi-acc-download -F fasta -m protein XP_021373098.1
blastp -db allprotein.fas -query XP_021373098.1.fa -outfmt 0 -max_hsps 1 -out opsin.blastp.typical.out


blastp -db allprotein.fas -query XP_021373098.1.fa -outfmt "6 sseqid pident length mismatch gapopen evalue bitscore pident stitle"  -max_hsps 1 -out opsin.blastp.detail.out

awk '{if ($6< 1e-35 )print $1 }' opsin.blastp.detail.out > opsin.blastp.detail.filtered.out
wc -l opsin.blastp.detail.filtered.out #48
                                               qseqid tophits    avPid avLength
3: Myessoensis.XP_021373098.1_rhodopsin_GQcoupledlike      48 36.29167 297.1042
   Myessoensis Lanatina Bbelcheri Aplanci Cintestinalis Avaga Dmelanogaster
3:           8        8         9       5             2     4             7
   Hsapiens Egranulosus minSp maxSp numZero unchar
3:        5           0     0     9       1  FALSE


grep -o -E "^[A-Z][a-z]+\." opsin.blastp.detail.filtered.out  | sort | uniq -c

   5 Aplanci.
   4 Avaga.
  10 Bbelcheri.
   1 Cintestinalis.
   7 Dmelanogaster.
   5 Hsapiens.
   8 Lanatina.
   8 Myessoensis.

grep -o -E "^[A-Z][a-z]+\." opsin.blastp.detail.filtered.out  | sort | uniq -c | wc -l 

8

 seqkit grep --pattern-file opsin.blastp.detail.filtered.out allprotein.fas > opsin.homologs.fas
muscle -in opsin.homologs.fas -out opsin.homologs.al.fas
 t_coffee -other_pg seq_reformat -in opsin.homologs.al.fas -output sim 

26.49% ID

alv -kli --majority opsin.homologs.al.fas | less -RS

iqtree -s opsin.homologs.al.fas  -m VT+F+R10
#started 2:04 - 2:11 (7 minutes)
37 dups, 45 losses (48 genes) score 100.5 (2.09points/gene)

# finding all groups of genes

#seqkit split2 --by-size 1 GCF_002113885.1_ASM211388v2.longest.renamed.faa


module load diamond/2.0.10
diamond makedb --in allprotein.fas --db allproteinDiamond

diamond blastp \
--db allproteinDiamond \
--query proteomes/GCF_002113885.1_ASM211388v2.longest.renamed.faa \
--very-sensitive \
--outfmt 6 qseqid sseqid pident length mismatch gapopen evalue bitscore pident stitle \
--evalue 1e-30 \
--max-hsps 1 \
--max-target-seqs 500 \
--threads 24  \
--out diamondFullMizuproteome.out 

module load R/4.1.1

library(data.table)
library(dplyr)
library(stringr)
library(matrixStats)

b1 <- fread("diamondFullMizuproteome.out")
colnames(b1) <- c("qseqid",	"sseqid",	"pident",	"length",	"mismatch",	"gapopen",	"evalue",	"bitscore",	"pident2",	"stitle")
b2 <- as.data.table(b1 %>% mutate(species = str_extract(stitle, "^[A-Z][a-z]+")))


myevalue <- 1e-30

b1e70 <- b2[evalue<myevalue, .(tophits = .N, avPid = mean(pident),avLength = mean(length), Myessoensis = sum(species=="Myessoensis"), 
Lanatina = sum(species=="Lanatina"), 
Bbelcheri = sum(species=="Bbelcheri"), 
Aplanci = sum(species=="Aplanci"), 
Cintestinalis = sum(species=="Cintestinalis"), 
Avaga = sum(species=="Avaga"), 
Dmelanogaster = sum(species=="Dmelanogaster"), 
Hsapiens = sum(species=="Hsapiens"), 
Egranulosus = sum(species=="Egranulosus")),
by = qseqid]

b1e70[, minSp := rowMins(as.matrix(.SD)), .SDcols = 5:13][] 
b1e70[, maxSp := rowMaxs(as.matrix(.SD)), .SDcols = 5:13][] 
b1e70[,  numZero := sum(.SD==0), .SDcols = 5:13,by=1:nrow(b1e70)][] 
b1e70[,  numZero := sum(.SD==0), .SDcols = 5:13,by=1:nrow(b1e70)][] 
b1e70[,  unchar := grepl("uncharacterized",.SD), .SDcols = 1,by=1:nrow(b1e70)][] 

###
b1opsin <- subset(b1e70,grepl("opsin",qseqid) & tophits <= 60 &  tophits >= 25 & numZero <= 1 & maxSp <= 15 & avLength > 200)

                                               qseqid tophits    avPid avLength
1:     Myessoensis.XP_021355729.1_rhodopsin_G0coupled      45 31.65333 302.6667
2:     Myessoensis.XP_021357501.1_rhodopsin_GQcoupled      59 35.38475 307.5424
3: Myessoensis.XP_021373098.1_rhodopsin_GQcoupledlike      48 36.29167 297.1042
   Myessoensis Lanatina Bbelcheri Aplanci Cintestinalis Avaga Dmelanogaster
1:          11        6        14       5             1     3             1
2:          11       10        14       5             1     6             7
3:           8        8         9       5             2     4             7
   Hsapiens Egranulosus minSp maxSp numZero unchar
1:        4           0     0    14       1  FALSE
2:        5           0     0    14       1  FALSE
3:        5           0     0     9       1  FALSE

###
b1e70F <- subset(b1e70,numZero == 0 & tophits <= 60 &  tophits >= 30 & maxSp <= 15 & minSp >= 2 & avLength > 200 & unchar==FALSE ) #

prohibitids <- ""
keepids <- ""

## now, go through hits. if any hit has self hits < evalue, remove it
for(dx in b1e70F$qseqid){
if(!(dx %in% prohibitids )){
keepids <- c(keepids,dx) 
prohibitids <- c(prohibitids,subset(b2,species=="Myessoensis" & evalue < myevalue & stitle == dx & qseqid != dx)$qseqid)
}
}
length(keepids) #134

save(b1e70F,file="b1e70Fnosacc.rda")
save(keepids,file="keepidsnosacc.rda")
write.table(keepids,file="idsb1e70nosacc.txt")

1e-39: 107
1e-30: 129


keep: Myessoensis.XP_021372857.1_solute_carrier_family_1
get rid of:
1: Myessoensis.XP_021372049.1_solute_carrier_family_1
2: Myessoensis.XP_021357580.1_solute_carrier_family_1
3: Myessoensis.XP_021354279.1_solute_carrier_family_1
4: Myessoensis.XP_021360761.1_solute_carrier_family_1

###### tags: `bio312`