HgcA.IMG.metagenomes
Robert W. Murdoch
December 11, 2018
Output summary
Metagenome search results
of 340449 proteins which were extracted from the IMG system, 12505 remained after hmmsearch (at 1e-20 threshold) and hhfilter
I am providing:
- fasta-file of all 340,000 proteins pulled out of IMG metagenome system
- pfam03599_all_mgenome.faa
- hmm results:
- hmm20result_e20.txt (report in tabular format)
- hmm20align_e20.stkhlm (alignment in stockholm format)
- gapped fasta formatted alignment:
- hmm20.align_e20.fasta
- hhfilter-ed alignment
- hgcA_IMG_mgenome.align.filtered.fasta
Genome search results
of 3574 proteins which were extracted from the IMG system, 551 remained after hmmsearch (at 1e-20 threshold) and hhfilter
I am providing:
- fasta-file of the 3574 proteins pulled out of IMG genome system
- pfam03599_all_genome.faa
- hmm results:
- hmm20result.g_e20.txt (report in tabular format)
- hmm20align.g_e20.stkhlm (alignment in stockholm format)
- gapped fasta formatted alignment:
- hmm20.align.g_e20.fasta
- hhfilter-ed alignment
- hgcA_IMG_genome.align.filtered.fasta
Ultimately, these can be traced back to particular scaffolds
Pipeline details:
1. All genes which contain pfam03599 from all publically available IMG metagenomes.
Genes were gathered from ~14,500 metagenomes in the IMG system.
Some broad information about the gene set
library(Biostrings)
pfam03599.all <- Biostrings::readAAStringSet("IMG_mgenome/seqs/pfam03599_all_mgenome.faa")
names <- names(pfam03599.all)
seqs <- paste(pfam03599.all)
length <- width(pfam03599.all)
pfam03599.all <- data.frame(names,length,seqs, stringsAsFactors = F)
kable(head(pfam03599.all))| names | length | seqs |
|---|---|---|
| 2001248046 CO dehydrogenase/acetyl-CoA synthase gamma subunit (corrinoid Fe-S protein) | 295 | MSQASCGCNSAKCQIVPVPVRTEPIAGTLTTEIAWRDVLGGWKARWGIGRMHYKISPGLYRVGQPSADSPVLVTANYKLTVDKVRKELGGMNVWLLILDTKGVNVWCAAGKGTFGIDELVRRVQAVNLKDIIKHRKLVLPQLGATGVAAHEVKRATGFSVHYGPVRARDIKAYLAAGMKATPAMRRVTFNWKERIVLTPVEIVGGMKWALVAAAVLAILDVIRHHRPTQHLAIDVLPFVAAVLTGGVLVPLCLPLLRSGAFALKGSVAGIVPIVALLLLLPTGGIPEPQEMRFSQ |
| 2001289001 | 73 | VLTAWAAGKFDAERIAKAVKSFNVADKIAAKRVVLPGAVAVLSGELEAELPGWEIRVGPREAVDLPKFMKQAL |
| 2001292656 CO dehydrogenase/acetyl-CoA synthase gamma subunit (corrinoid Fe-S protein) | 255 | VARAARAPLAVRAASLDQLATLTEALKKAGVEDLVIDPLAREPGDSLARLTQIRRLAVKKNFRALGYPVITFPGEATGSDSALLAAQHIAKYAGFIVLERFDPADLYPLLVLRQNLYTDPQKPIQVQPGLYEINAPTADAPVMVTTNFSITYFAVSNEVEGSGRPGWLVVTDAEGMSVLTAWAAGKFDADRIAKAVKTYGVAERVNHRRLILPGQVAVLSGEVEESLPGWEIMVGPKEAVDLPGYLKLVWGQPTP |
| 2001322166 CO dehydrogenase/acetyl-CoA synthase delta subunit (corrinoid Fe-S protein) | 191 | MPIEIPKDKWPGSIKAITIGATAAEGGTRTRTVTVGGQTTMPYLQYEAPTPNKPVIAIEIKSRKPEDWSALLAETWGAAMTDPAQWAQKAEEAGADLIVLTLTVEDSVEDAVRVVKSVLAATGLPIAVWGPGQAEKDNELLVPVSEATKGEKILLGICEDKNYRTIVATAMANGHLVQSRAPMDGTIEHNP |
| orf47_0114_11 CO dehydrogenase/acetyl-CoA synthase delta subun it | 525 | LAFEFFKESYSGNIKEITLGKGDSAVTVGGEACYPFYQFEGAMPNKPKIAMEIWDMVPDEWPEAALSPFKDVVSDPAAWAKKCVEEYGAEMIVLQLKSTDPNGENASADDAVASVKKVLGAVDVPLIIWGTANIEKDEEVLKQIAEECQGENLILGPVEDKNHKGIGAAAMGFGHSIIASSPIDVNLAKQCNILLENLGMPLDRMIIDPTTGGLGYGLEYSYSVMERIRMAALAQGDDKLQLPLINNLGNEVWKCKEAKQSVEEAPVLGDPEKRPILMEAVGAVAYLMSGSNILIMRHPESIRLVKAYIDLVAEGGLASDVAPIAKQLDDVEIDFLALSPDPDLTIPEEKEKKPAPAKKAKAPAPAKEAAPAAAKAEPAQKAAPAAKPAAPAAAPEVDAAADEKAKAEAAAKAKAEAEAKAKLDAAAQAKADEEARIKADAEAKAKAEADAKAKADEAARQAAEEKAKVEAELDSVREKRAKEREALAAKRQSAEEDEVSKTPAEVQKSMTDKMVENLDWVHRRI |
| orf51_0114_11 CO dehydrogenase/acetyl-CoA synthase delta subun it:Putative Fe-S cluster | 446 | MALTGIQIFKLLPKTNCKECGVPTCLAFAMNLASGKAELDSCPYVSDEAREQLSEASAPPIRPVAIGKGVRAATSGGETVQYRHEKTFFNPTLLAATVGSDIAEADLETKLKAWNAFQFERVGLNLRPELIVLKDAGGNGDAFAKTAKMIAETSEFNLALMTEDVGVMKAGIEACGFKRPLMYAATEANVDDYCALAKENDLPLAVKADSVEGLVPLTEKLTEAGLKDLVIDPGSREIKQAHEDQIALRRAALKDGNRALGFPTITFPCEMASNLDMETLIAGTFVAKYGGIAVLSDFAGESLFPLLLERLNIFTDPQRPMTVTEGIYEVNNPDENSPVMVTTNFALTYFIVAGEIEGSRVPAWLLIKDSEGLSVMTAWAAGKFAGDDVGMFVKKCGIADKIKHQELIIPGYAAAIAGDVEEELPGWTITVGPREAAHIPGFLKAK |
## ## There were 340449 proteins retrieved from publically available IMG metagenomes.
##
## The mean protein length was 194.467224165734
## The median length was 149boxplot(pfam03599.all$length, horizontal = TRUE, main = "Protein Length Distribution",
xlab = "protein length (amino acids)")
2. Apply HgcA hmm model generated from the uniprot20 database to the gene set.
Code used to query the model:
$ hmmsearch -A hmm20align.stkhlm –tblout hmm20result.txt –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa
The search was done using default settings. What should be the threshold for inclusion?
hmmsearch results
hmm.result <- read.table("IMG_mgenome/seqs/hmm20result.txt", fill=TRUE, stringsAsFactors = F)
kable(head(hmm.result))| V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | V11 | V12 | V13 | V14 | V15 | V16 | V17 | V18 | V19 | V20 | V21 | V22 | V23 | V24 | V25 | V26 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Ga0066906_100185416 | - | hgca.uniprot20 | - | 7.4e-146 | 489.5 | 0.8 | 8.5e-146 | 489.3 | 0.8 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
| Ga0153798_100215862 | - | hgca.uniprot20 | - | 7.4e-146 | 489.5 | 0.8 | 8.5e-146 | 489.3 | 0.8 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
| Ga0134093_10030757 | - | hgca.uniprot20 | - | 2.7e-145 | 487.6 | 2.0 | 3.2e-145 | 487.4 | 2.0 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
| Ga0134094_10027597 | - | hgca.uniprot20 | - | 2.7e-145 | 487.6 | 2.0 | 3.2e-145 | 487.4 | 2.0 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
| Ga0134095_10158963 | - | hgca.uniprot20 | - | 2.7e-145 | 487.6 | 2.0 | 3.2e-145 | 487.4 | 2.0 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
| Ga0134097_100290211 | - | hgca.uniprot20 | - | 2.7e-145 | 487.6 | 2.0 | 3.2e-145 | 487.4 | 2.0 | 1.0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | CO | dehydrogenase/acetyl-CoA | synthase | gamma | subunit | (corrinoid | Fe-S | protein) |
## under default conditions, 141025 proteins were in the reportThe problem here is that the default reporting setting is extremely lenient. Let’s run the search again with 0.01, the default inclusion threshold for the MSA generation:
Drop reporting threshold from e=10 to e=0.01
$ hmmsearch -A hmm20align_e0.01.stkhlm –tblout hmm20result_e0.01.txt -E 0.01 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa
hmm.result0.01 <- read.table("IMG_mgenome/seqs/hmm20result_e0.01.txt", fill=TRUE, stringsAsFactors = F)## under default conditions, 134780 proteins were in the reportThis changed almost nothing at all.
Drop reporting threshold to e=0.00001 (e^-5)
$ hmmsearch -A hmm20align_e5.stkhlm –tblout hmm20result_e5.txt -E 0.00001 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa
hmm.result_e5 <- read.table("IMG_mgenome/seqs/hmm20result_e5.txt", fill=TRUE, stringsAsFactors = F)## under default conditions, 121488 proteins were in the reportAgain, almost no difference.
complete a cutoff threshold series down to e^-30
hmmsearch -A hmm20align_e10.stkhlm –tblout hmm20result_e10.txt -E 0.0000000001 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa hmmsearch -A hmm20align_e20.stkhlm –tblout hmm20result_e20.txt -E 0.00000000000000000001 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa hmmsearch -A hmm20align_e30.stkhlm –tblout hmm20result_e30.txt -E 0.000000000000000000000000000001 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa
hmm.result_e10 <- read.table("IMG_mgenome/seqs/hmm20result_e10.txt", fill=TRUE, stringsAsFactors = F)
hmm.result_e20 <- read.table("IMG_mgenome/seqs/hmm20result_e20.txt", fill=TRUE, stringsAsFactors = F)
hmm.result_e30 <- read.table("IMG_mgenome/seqs/hmm20result_e30.txt", fill=TRUE, stringsAsFactors = F)## e.series homologs
## 1 1e+01 141025
## 2 1e-02 134780
## 3 1e-05 121488
## 4 1e-10 99264
## 5 1e-20 53346
## 6 1e-30 387813. Filtering the original MSA using hhfilter
convert the stockholm format MSA to gapped fasta using https://biocs-blog.blogspot.com/2010/08/convert-stockholm-sequence-format-to.html?m=1
install hh-suite from https://github.com/soedinglab/hh-suite
filter the alignment using default parameters via:
$ hhfilter -i hmm20.align.fasta -o hmm20.align.filtered.fasta
Received a warning message: WARNING: Maximum number 65535 of sequences exceeded in file hmm20.align.fasta
The gapped alignment has 132267 sequences (why is this less than the e=0.01 report?). However, these are listed in the file in order of e-value, same as the order that appears in the report. Thus, we can use a stricter alignment construction setting.
Use the hmm stockholm alignment with 1e-20 cutoff, convert to fasta, and filter:
$ hmmsearch -A hmm20align_e20.stkhlm –tblout hmm20result_e20.txt -E 0.00000000000000000001 –notextw –cpu 6 ../hgca.uniprot20.hmm 58401.faa $ perl stock_2_fasta.pl hmm20align_e20.stkhlm -g > hmm20.align_e20.fasta $ hhfilter -i hmm20.align_e20.fasta -o hmm20.align_e20.filtered.fasta
library(Biostrings)
hhfiltered <- Biostrings::readAAStringSet("IMG_mgenome/seqs/hmm20.align_e20.filtered.fasta")
names <- names(hhfiltered)
seqs <- paste(hhfiltered)
length <- width(hhfiltered)
hhfiltered <- data.frame(names,length,seqs)
kable(head(hhfiltered))| names | length | seqs |
|---|---|---|
| Ga0066906_100185416/7-334 | 332 | –SREIKIMQTSDQLTFKDKLGAWKARWGINRMNYKVEPGLYRLGNPDGNSPVLVTANYKMTFDSLRKELSGVDAWILVLDTKGINVWCAAGKGTFGTRELINRLAIVQLEKVINHRTLILPQLGAPGISAHEVTKNTGFKVVYGPVRAKDIKEFIRSGMQATPEMRTVKFSAYDRLVLTPIEFVGTIKPSLMIFGVLFLLNLMGLG–PFGLVDFYGYVGAVIIGCILTPVLLPLIPGRAFAWKGWLLGLLWAVLVNWLngWpsnPQYGLLQALGYVLVFPAISAYLAMNFTGSSTYTSLSGVLKEMKIAIPAIIISIVLGCLLLLINSFI |
| Ga0134097_100290211/8-334 | 332 | —REIKIMQTSPQLTFRDKLGAWKVRWGINRMNYKVEPGLYRLGNPGENSPVLVTANYKMTFDSLRKELSELDAWILVLDTKGINVWCAAGAGTFGTKELINRMAITQLYKVVTHRTLILPQLGAPGVSAHEVKKYTGFKVVYGPVRAKDLKRFIRSGMKATPEMRTVKFNVYDRLVLTPVELVSTVKPSLMIFGVLFLLNLIGLG–PFGLVDFYGYVGSVIIGCVLTPVLLPLIPGRAFAWKGWLLGFIWALLVNWLNgwpfnPQYSLLQALGYVLVLPAISAYLALNFTGSSTYTSLSGVLKEMRIAIPAIIVSIVIGCLLILVNNFI |
| Ga0134099_10091753/2-322 | 332 | ———QTSSQLTFQDKLGAWKARWGINRMNYKVEPGLYSLGNPDENSPVLVTANYKLTFDSLRKELSGLDAWILVLDTKGINVWCAAGKGTFGTTELVNRLLVVQLEKFVAHRTLILPQLGAPGVSAHEVTKLTGFKVVYGPVRAKDIKEFIQSGMKATPEMRTVKFSAYDRLVLTPIELVGTIKPSLMIFGVLFLLNLMGLG–PFGLVDFYGYMGAVIIGCVLTPVLLPLIPGRAFAWKGWLLGFVWAVLVNVLNgwpddPQYSLLRAIGYLLVFPAVSAYLALNFTGSSTYTSLSGVLKEMRIAIPVIIISIVIGCLLILINRFI |
| Incfw_10085297/71-398 | 332 | -KTDAGLVPQVSTKLVADDILGAWKVRWGINRMNYKINPGIYAVGTPDKDSTVLVTANYKLTFDALRRELQGLNAWILVLDTKGINVWCAAGKGTFGTKELINRISKTKLSQIVTHKNLILPQLGGPGVSAHEITKHTGFKVVYGPVRARDIKEFIKSGMKATDEMRAVKFNTYDRLVLTPMELVGTFKISLIIFGILFLLNLVAL–IPFGIVDFYAYAGALFTGCVLTPVLLPWIPGRAFAWKGWLLGFIWAVAVNILngWpetPSYSLIRAIGYLLILPSVSSFYAMNFTGSSTYTSLSGVLKEMKIAVPVIIISISTGIILVLLNIF- |
| Ga0134096_10031656/20-342 | 329 | –TPKGSVPVISTNLHFKDILGAWKVRWGIGRMNYKISPGLYAIGKPDHTSPVLVSANYKLTFDVLRKELSGLDCWILILDTKGINVWCAAGKGTFGTDELVNRISKTGLSEIVMHRKLVLPQLGAPGVSAHEVTKQTGFSVVYGPVRAKDIKAFLDSGFKATEEMRTVKFTMWDRLVLTPVELVAAAKFSLMVFGVMFLLNLFAAS–PFGLADFIAYTGAVVTGAVITPMLLPMIPGRAFAWKGWLLGLLWTLgfIWYSGWLVSELLLAMGYLLVLPSLSAFLAMNFTGSSTYTSFSGVIKEMKIAVPLIALSTIAGIVLVLIQS– |
| Ga0153798_100569682/35-360 | 330 | -HTPKGSVPVVSTNLHFKDILGAWKVRWGIGRMNYKINPGLYAIGKPDHTSPVLVSANYKLTFDALRKELSDLDCWILILDTKGINVWCSAGKGTFGTDELVNRISKTGLSEIVSHRKLVLPQLGAPGVSAHEVTKQTGFSVVYGPVRAKDIKAFFDSGFKATEEMRTVKFTMWDRLVLTPVELVAAAKFSLMVFGIMFLLNLFAAR–SFGLADFIAYAGTVVAGTVITPLLLPVIPGRAFAWKGWLLGLLWALGFIWYSgwfVPESLLLAIGYLLVLPSFSAYLTMNFTGSSTYTSFSGVIKEMKIAVPLIAISLIVGIVLLLINSI- |
## # There are 12505 proteins in the final filtered set, with an average length of 330.0934boxplot(hhfiltered$length, horizontal = TRUE, main = "Protein Length Distribution",
xlab = "protein length (amino acids)")
4. Same pipeline applied to genomes only
library(Biostrings)
pfam03599.genomes <- Biostrings::readAAStringSet("IMG_genome/seqs/pfam03599_all_genome.faa")
names <- names(pfam03599.genomes)
seqs <- paste(pfam03599.genomes)
length <- width(pfam03599.genomes)
pfam03599.genomes <- data.frame(names,length,seqs, stringsAsFactors = F)
kable(head(pfam03599.genomes))| names | length | seqs |
|---|---|---|
| 637778520 acetyl-CoA decarbonylase/synthase complex subunit gamma | 323 | MPIVSTEFRFADRLGAWKARWGIGRMSYLVPPGLYAIGSPGPEDPVVVTANYKMSYDLVRRELSGRNVWFLVLETFGINVWCAAGKGTFGTDELVRRVGMVKLATVVKHRALILPILGAPGVAAHQVARRSGFRVHYATIRAADLPEYLDNGMVTTPEMRQLTFTTRERLALTPIELVTTGKWVVPFVALLVVLSIVPGGLLSSPTLVSSLVLFVGTVLAGAVAAPVLLPWLPGKSFAVKGAVAGCCWAVLFIGVGRGATGWLDAMALFLTVPAVTAFLTLNFTGSTPFTSRSGVKKEMRLGLPMMAVSLLAGISLWIAARFF |
| 637861274 carbon monoxide dehydrogenase corrinoid/iron-sulfur protein gamma subunit | 396 | MMKDSEKPGKSDMSEQRKPVMIFPMIQPGFTAGRMGREVPLNLPAITQSFVVDTLSTAVGMVPVVDSRLTRQDHVGSLKARFGIGRMNYRVDPGLYALGRPDESAPVLVSANYKLSFDCLREALPGRTAWILVLDTNGINVWCAAGKGAFGTEEVVRRLRSSGLEKLVKHRRLILPQLGASGVAAHEVRKQTGFAVSYGPVEARDLPAYLDAGMKATEAMRRKNFPFRDRLVLIPMECLPALKPLFLLTLFFLIAGGIVQGPRWTVEPTLNSGIFMTFALLGAFFAGTIFTPAFLPWLPGRAFSAKGLCPAWIVVLILLESRLTFPLSLPRVLESFSVLFLVPAAASYWAMKFTGSTPITSLSGVRKEMRRALPLQIAGALIGLGCWAAAKLLRVL |
| 637958648 acetyl-CoA decarbonylase/synthase gamma subunit (EC 2.1.1.-) | 468 | MKINSPLEAYKFLPGTNCGECGETSCMAFASHLIDRSLKATDCTPLVNESKYKKKYDELVTLLAPEIREVVIGVGDNAISIGGDDVLHRHKLTFFNKTALAYDVWDTMEETDLVERVNKIQDFKKFYVGDYLTLDMIAVRSVSNDPEKFAATVKKVMETTEFPMILCSFEPAVLRAGLEVAAEKRPLLYAANKDNWQEVAALAQEFNVPVTVFAPNDLDLLKSMAKTFADNAIEDVVLDPGTFPTGKGLKTTLNNFLKIRRASINGDRDIAFPIMAVPFTAWMAHDDEVSASYWETVVASIFTIKYGDIMIMHSIEPYALMPQLHIRDTMYTDPRKPVTVDPAVYEIGSPTADSPLLVTTNFALTYYTVESDLSSNGVDCYLTAIDTDGIGVEAAVAGGQLTAAKIKKGLDDAGFDMEKLTHKVIVLPGLAARLQGDVEDETGAGVMIGPADSGRLPGWMEQNWPPQK |
| 637958647 acetyl-CoA decarbonylase/synthase delta subunit | 439 | MSNKMKLSGLTDILKDLDVQSLEGVTIEGDIELDISGGGGLNPAIAYALGHEAAQISMHVANIARMLGYPVDQLFAASMGGVPQQAPVAGASRIPELLPAKFDVSKMSEWATPIQEVTLGATSADGGSRKSTVTLGGENALPYYFDAEMPHRNYVTMDVFDMPISMAKSVKENYGDVINDPAEWAKKVVKDFNADMVTIHLISTDPLINDTPAREAAKVVEDVLQAVKVPIVIGGSGNPKKDPEVLERAAEVAEGERVLLASASLNLDYERVAKAATDHGHAVLSWTQLEINAQKELNRKLMKQCNVPRDSIVMDPTTAALGYGLDYAYTNMERIRLAGLMGDDELTFPMSSGTTNAWGARESWMKESPLKQDSDWGPREYRGPIWEIITGLTLSLAGNDMFMMMHPTSVQVLKEITQTLYGSIEAEEIDITKWIGAEV |
| 637896512 hypothetical protein | 342 | METSCCCSCSEPDTACSYSIPQTDSNITIANRIDHILARWGYKRSGHVVRPGLYRLGNPDADSPVFASANYTISFDTLRLNLKGFDAYILVLDTKGINVWCAAGKGTFGTDELVKRIESTGLSHIVSHRNIIVPQLGAPGISAHMVKQRSGFTVEYGPVRAADLPEYMRTGKATSEMRKVSFPIQDRLVLIPIELVHVAIPTIIISIILWLLAGPAAGVAAIAAVITGTVLFPILLPYLPTQDFSMKGFILGGVCAVPFSAMYAQSADLPEWALAAGAIVPLLIIPAVVAYLGLNFTGSTPYTSRTGVKKEIFRYVPIMAVMTGIGSVIGIVLSIFRLAGVI |
| 637896324 acetyl-CoA decarbonylase/synthase delta subunit | 416 | MTDQEKPDISNLIRLLGPGLGELLSGRQIELEHVELELGELELWIPVGGISQGFPAPVAAAAKQMTPPVKRLPRPDDLISEKFLMEPDVFPAQIREITLGATRAEGGTRGSTITVGGSTSIPFTHPQSPPPHPPVISLDVFDTRVPMPKALKNHLTEVIDDPAAWAKMNVEQFGADMVTVHLLSTDPLIQNKSPKEASKTVEEVLQAVDVPLIIGGCGDPKKDAETFTEIAAMAEGERLLLNSVTSDMADAKTLEPVCKAADTHGHCLLGFTGLDLNSAKELNRRIYQYFPPERLLMDLTTVALGYGLEYSFSIHERARMAALMGDPELQHPTISACTNAWSAREAWMDLGPAFGRNDLRGPLWETIGGITLLLAGVDLFLMMHPLAVSTLRDVAGKLMHPAEIKPDMADWAALKI |
## ## There were 3574 proteins retrieved from publically available IMG genomes.
##
## The mean protein length was 388.393116955792
## The median length was 402boxplot(pfam03599.genomes$length, horizontal = TRUE, main = "Protein Length Distribution",
xlab = "protein length (amino acids)")
hmmsearch -A hmm20align.g_e20.stkhlm –tblout hmm20result.g_e20.txt -E 0.00000000000000000001 –notextw –cpu 6 hgca.uniprot20.hmm pfam03599_all_genome.faa perl /home/rmurdoch/Tools/stock_2_fasta.pl hmm20align.g_e20.stkhlm -g > hmm20.align.g_e20.fasta hhfilter -i hmm20.align.g_e20.fasta -o hgcA_IMG_genome.align.filtered.fasta
library(Biostrings)
hhfiltered.g <- Biostrings::readAAStringSet("IMG_genome/seqs/hgcA_IMG_genome.align.filtered.fasta")
names <- names(hhfiltered.g)
seqs <- paste(hhfiltered.g)
length <- width(hhfiltered.g)
hhfiltered.g <- data.frame(names,length,seqs)
kable(head(hhfiltered.g))| names | length | seqs |
|---|---|---|
| 2507427874/8-334 | 332 | —REIKIMQTSTQLTFHDRLGSWKARWGINRMNYKVEAGLYRVGDPNENSPVLVTANYKMSFDSLRKELSGLDTWILVLDTKGINVWCAAGKGTFGTTELLNRLAIVQLDRVISHRTLILPQLGAPGVSAHEITKYSGFKVIYGPVRAKDLKEFIQAGMKATPEMRKVRFSAYDRLVLTPIELVGTLKASLMIFGILFLLNLMGLG–PFGLVDFYGYMGAVIIGCVLTPVLLPWIPGRAFAWKGWLLGLVWTVFVNATngWpydPQYSMLQALGYVLVFPAVSAYLAMNFTGSSTYTSFSGVLKEMRIAIPAIIVSIFIGCVLILVNSFI |
| 2510774674/71-398 | 332 | -KTPRGKVPIVSTQLAFSDRLGGWKARWGINRMNYKINPGLYAVGIPDNTSPVLVTANYKLTFDALRKELADINAWILVLDTRGINVWCAAGKGTFGTREIVTRINKTRLHEIVSHRTLILPQLGAPGVSAHEVSTKIGFKVVYGPVRASDIKEFITLGMKATPEMRTVKFNTYDRLVLTPMELVNTFKICLMVLGLLFILNLIAVR–SFGVIDLYAYMGAILAGTVLTPVLLPWIPGRAFAWKGWFLGFIWAITVNILngwpfSPSYSLIRLLGYLLILPSISSFYAMNFTGSSTYTSLSGVLKEMKIAIPIIIITISSGVLLILIDSF- |
| 2571048919/52-376 | 332 | —–DSITKVSTKLTFKDTLGGWKARWGINRMNYKVNPGIYAVGDPDKSSPVLVTANYKMTFDALRKELSGLSAWIMVLDTKGINVWCAAGKGTFGTKELITRIGKIRLNDIVAHRTLILPQLGAPGISAHEVTKSTGFKIVYGPVRATDIKEFIRSGMKATGEMRTVKFSTYDRLVLTPMELVPTFKVSMIVFGILFLLNLIAAK–PFGMVDFYGYAGAVITGCVLAPVLLPWIPGKAFAWKGWLLGFLWAIAVNIIngwpqLPSYSLVRALGYMLILPAVSSFYAMNFTGSSTCTSFSGVLKEMKMAIPIIIAMIALGAVMILLNSFI |
| 2508509747/5-334 | 334 | –SRELNITPTSAQLTFSDILGTWKARWGIGRMKYMVEPGLYSLGKPDSNSPVLVSANYKMSFDKLRKELKGIDAWILVLDTKGINVWCAAGKGTFGTQELLNRMAIVQLEKVVSHRTVIVPQLGAPGISAHVVAKFSGFKVVYGPVRAEDIKQFIKSGMKATQEMRTVRFGIYDRLVLTPVELVGTFKISLMIFGVLFLLNLLGLG–PFGLVDFYAYMGAVVAGCVLTPVLLPWIPGRAFAWKGWVLGFIWAIIVNILndwngtgASEYSSLRALGYLFLLPSISAFLAMNFTGSSTFTSFSGVLKEMRKAVPAIIISTVLGMVLILVNSFM |
| 2507227236/7-333 | 332 | –KKEIKILQTSSQLTFADTLGAWKVRWGIDRMNYKVEPGLYRVGNPDGSSPVLVTANYKLTFDGLRKKLTGLDTWILALDTKGVNVWCAAGKGTFGTTELINRLAIAQLDKVVAHRTLILPQLGAPGVSAHEVTKYTGFRVIYGPVKAEDIKEFIQAGMKATPEMRTVQFGAYDRLILTPIELVGTLKKSLMIFGVLFLLNLMGLG–PFGWVDFYGYIGAVIIGCVLTPALLPWIPGRAFAWKGWLLGLVWAILVNGLngWpyePRYSLLQALGYILIFPAVSGYLAMNFTGSSTYTSFSGVLKEMRIAIPTIIISIIVGCVLLLAGSF- |
| 2509272264/8-334 | 332 | —REIKITPTSSQLTAQDIRGTWKARWGIGRMNYKVEPGLYRLGNPDPNAPVLVTANYKMTFDALRKELTDLDAWILVLDTKGVNVWCAAGKGTFGTTELINRLAVVQLEKVVSHRLLILPQLGAPGVSAHEVLKFSGFKVLYGPVRAKDIKEFIRSGMKATGDMRKVKFNTYDRIVLTPVELVITFKKSFVIFGVLFLLNLMGIG–PFGLVDFYSYLGAVIIGCVLTPVLLPWIPGRAFAWKGWLLGFIWAVLVNFLngwpaAPEYSLLRALGYMLILPSVSAYLAMNFTGSSTYTSFSGVLKEMKIAIPGIIISLVLGAVLILINSFI |
## # There are 551 proteins in the final filtered set, with an average length of 331.1434boxplot(hhfiltered.g$length, horizontal = TRUE, main = "Protein Length Distribution",
xlab = "protein length (amino acids)")