# 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`