text:
heatmap of specific sequence motifs in aligned fasta files
I have a collection of fasta files, each containing three aligned sequences. I am interested in understanding the distribution of a specific sequence motif PGP, RGP and KGP in all the alignments. I was wondering if there is possible way in R to create a heatmap highlighting the location of these motifs.
I have tried the code in this post to create a heatmap but it is not working as expected. One issue is that ggmsa does not import the fasta file. The fasta file has to be copied into the installation directory for it to work. I have close to 50 aligned files, so not very efficient. The other major problem is it creates a heatmap compared to a reference sequence and not the motifs. Can't figure how to change this.
Any help is much appreciated. Two fasta files containing aligned sequence are as follows,
>chain A CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain B CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain C CO1A2_HUMAN/1-1366 Collagen alpha-2(I) chain
----------------------------------------------------------------------------------------MLSFVDTRTLLLLAVTLCLATCQSLQEETVRKGPAGDRGPRGERGPPGPPGRDGEDGPTGPPGPPGPPGPPGLGGNFAAQYDGKGVGLGPGPMGLMGPRGPPGAAGAPGPQGFQGPAGEPGEPGQTGPAGARGPAGPPGKAGEDGHPGKPGRPGERGVVGPQGARGFPGTPGLPGFKGIRGHNGLDGLKGQPGAPGVKGEPGAPGENGTPGQTGARGLPGERGRVGAPGPAGARGSDGSVGPVGPAGPIGSAGPPGFPGAPGPKGEIGAVGNAGPAGPAGPRGEVGLPGLSGPVGPPGNPGANGLTGAKGAAGLPGVAGAPGLPGPRGIPGPVGAAGATGARGLVGEPGPAGSKGESGNKGEPGSAGPQGPPGPSGEEGKRGPNGEAGSAGPPGPPGLRGSPGSRGLPGADGRAGVMGPPGSRGASGPAGVRGPNGDAGRPGEPGLMGPRGLPGSPGNIGPAGKEGPVGLPGIDGRPGPIGPAGARGEPGNIGFPGPKGPTGDPGKNGDKGHAGLAGARGAPGPDGNNGAQGPPGPQGVQGGKGEQGPPGPPGFQGLPGPSGPAGEVGKPGERGLHGEFGLPGPAGPRGERGPPGESGAAGPTGPIGSRGPSGPPGPDGNKGEPGVVGAVGTAGPSGPSGLPGERGAAGIPGGKGEKGEPGLRGEIGNPGRDGARGAPGAVGAPGPAGATGDRGEAGAAGPAGPAGPRGSPGERGEVGPAGPNGFAGPAGAAGQPGAKGERGAKGPKGENGVVGPTGPVGAAGPAGPNGPPGPAGSRGDGGPPGMTGFPGAAGRTGPPGPSGISGPPGPPGPAGKEGLRGPRGDQGPVGRTGEVGAVGPPGFAGEKGPSGEAGTAGPPGTPGPQGLLGAPGILGLPGSRGERGLPGVAGAVGEPGPLGIAGPPGARGPPGAVGSPGVNGAPGEAGRDGNPGNDGPPGRDGQPGHKGERGYPGNIGPVGAAGAPGPHGPVGPAGKHGNRGETGPSGPVGPAGAVGPRGPSGPQGIRGDKGEPGEKGPRGLPGLKGHNGLQGLPGIAGHHGDQGAPGSVGPAGPRGPAGPSGPAGKDGRTGHPGTVGPAGIRGPQGHQGPAGPPGPPGPPGPPGVSGGGYDFGYDGDFYRADQPRSAPSLRPKDYEVDATLKSLNNQIETLLTPEGSRKNPARTCRDLRLSHPEWSSGYYWIDPNQGCTMDAIKVYCDFSTGETCIRAQPENIPAKNWYRSSKDKKHVWLGETINAGSQFEYNVEGVTSKEMATQLAFMRLLANYASQNITYHCKNSIAYMDEETGNLKKAVILQGSNDVELVAEGNSRFTYTVLVDGCSKKTNEWGKTIIEYKTNKPSRLPFLDIAPLDIGGADQEFFVDIGPVCFK

and
>chain A CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain B CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain C CO4A2_HUMAN/1-1712 Collagen alpha-2(IV) chain
MGRDQRAVAGPALRRWLLLGTVTVGFLAQSVLAGVKKFDVPCGGRDCSGGCQCYPEKGGRGQPGPVGPQGYNGPPGLQGFPGLQGRKGDKGERGAPGVTGPKGDVGARGVSGFPGADGIPGHPGQGGPRGRPGYDGCNGTQGDSGPQGPPGSEGFTGPPGPQGPKGQKGEP-YALPKEERDRYRGEPGEPGLVGFQGPPGRPGHVGQMGPVGAPGRPGPPGPPGPKGQQGNRGLGFYGVKGEKGDVGQPGPNGIPSDTLHPIIAPTGVTFHPDQYKGEKGSEGEPGIRGISLKGEEGIMGFPGLRGYPGLSGEKGSPGQKGSRGLDGYQGPDGPRGPKGEAGDPGPPGLP--AYSPHPSLAKGARGDPGFPGAQGEPGSQGEPGDPGLPGPPGLSIGDGDQRRGLPGEMGPKGFIGDPGIPALYGGPPGPDGKRGPPGPPGLPGPPGPDGFL-FGLKGAKGRAGFPGLPGSPGARGPKGWKGDAGECRCTEGDEAIKGLPGLPGPKGFAGINGEPGRKGDRGDPGQHGLPGFPGLKGVPGNIGAPGPKGAKGDS-RTITTKGERGQPGVPGVPGMKGDDGSPGRDGLDGFPGLPGPPGD-GIKGPPGDPGYPGIPGTKGTPGEMGPPGLGLPGLKGQRGFPGDAGLPGPPGFLGPPGPAGTPGQIDCDTDVKRAVGGDRQEAIQPGCIGGPKGLPGLPGPPGPTGAKGLRGIPGFAGADGGPGPRGLPGDAGREGFPGPPGFIGPRGSKGAVGLPGPDGSPGPIGLPGPDGPPGERGLPGEVLGAQPGPRGDAGVPGQPGLKGLPGDRGPPGFRGSQGMPGMPGLKGQPGLPGPSGQPGLYGPPGLHGFPGAPGQEGPLGLPGIPGREGLPGDRGDPGDTGAPGPVGMKGLSGDRGDAGFTGEQGHPGSPGFKGIDGMPGTPGLKGDRGSPGMDGFQGMPGLKGRPGFPGSKGEAGFFGIPGLKGLAGEPGFKGSRGDPGPPGPP-PVILPGMKDIKGEKGDEGPMGLKGYLGAKGIQGMPGIPGLSGIPGLPGRPGHIKGVKGDIGVPGIPGLPGFPGVAGPPGITGFPGFIGSRGDKGAPGRAGLYGEIGATGDFGDIGDT-INLPGRPGLKGERGTTGIPGLKGFFGEKGTEGDIGFPGITGVTGVQGPPGLKGQTGFPGLTGPPGSQGELGRIGLPGGKGDDGWPGAPGLPGFPGLRGIRGLHGLPGTKGFPGSPGSDIHGDPGFPGPPGERGDPGEANTLPGPVGVPGQKGDQGAPGERGPPGSPGLQGFPGITPPSNISGAPGDKGAPGIFGLKGYRGPPGPPGSAALPGSKGDTGNPGAPGTPGTKGWAGDSGPQGRPGVFGLPGEKGPRGEQGFMGNTGPTGAVGDRGPKGPKGDPGFPGAPGTVGAPGIAGIPQKIAVQPGTVGPQGRRGPPGAPGEMGPQGPPGEPGFRGAPGKAGPQGRGGVSAVPGFRGDEGPIGHQGPIGQEGAPGRPGSPGLPGMPGR-SVSIGYLLVKHSQTDQEPMCPVGMNKLWSGYSLLYFEGQEKAHNQDLGLAGSCLARFSTMPFLYCNPGDVCYYASRNDKSYWLSTTAPLP--MMPVAEDEIKPYISRCSVCEAPAIAIAVHSQDVSIPHCPAGWRSLWIGYSFLMHTAAGDEGGGQSLVSPGSCLEDFRATPFIECNGGRGTCHYYANKYSFWLTTIPEQSFQGSPSADTLKAGLIRTHISRCQVCMKNL

When I want an MSA visualization in R that I can't find pre-made, my approach is to put the text into a data frame of position/character combos and use generic plotting tools like pheatmap or ggplot2's geom_tile.  Looks like that'd get tricky here because of the gaps and multi-AA motifs though.  Maybe it'd work with a mapping of gapped/ungapped locations of all AA trios instead of individual characters?  And for simple FASTA reading maybe check out this minimalistic package from a colleague of mine.  It just reads the text into data frames.

Here's a solution based on Jesse's comment. The steps are:

read in MSA with Biostrings::readAAStringSet()
right-fill sequences so they are all the same length (assumes left alignment)
search for motif substrings in each amino acid sequence with gregexpr(). If a motif in the motif set is found, record position of center of motif.
bind the results from each sequence to a single table

find_motifs <- function(msa, motifs) {

  require(tibble)
  require(Biostrings)

  # read in fasta as Biostrings object
  seqs <- Biostrings::readAAStringSet(filepath = msa, format = "fasta")
  maxseqlen <- max(width(seqs))

  # initiate results tibble
  results <- tibble::tibble(position = seq(1:maxseqlen))

  # readAAMultipleAlignment requires equal length of sequences
  # instead, use readAAStringSet and fill "-" on right
  # assumes left-aligned MSA
  for (i in 1:length(seqs)) {

    seqlen <- nchar(seqs[i])

    if (seqlen < maxseqlen) {
      fill <- maxseqlen - seqlen
      seqs[[i]] <- paste0(as.character(seqs[[i]]),
                          paste(rep("-", fill), collapse = ""))
    }

    # initiate vector of zeros to store hit positions
    positions <- rep(0, maxseqlen)

    lapply(motifs, function(motif) {

      # factor to add to gregexpr output to center hit
      # for motifs with an even number of characters, use left center
      shift <- ceiling(nchar(motif)/2)

      hits <- unlist(gregexpr(motif, seqs[[i]])) + shift

      lapply(hits, function(hit) {
        positions[hit] <<- 1
      }) # end apply over hits

    }) # end apply over motifs

    results <- dplyr::bind_cols(results, positions)

  } # end loop over sequences

  names(results) <- c("position", paste0("seq", seq(1:length(seqs))))
  results

} # end function

Running this function on the provided fastas, we get a tables where positions in sequences matching motifs are recorded as 1.
setwd("/path/to/files")
ngp <- c("PGP","RGP","KGP")
msa1 <- find_motifs(msa = "msa1.fasta",
                    motifs = ngp)
msa2 <- find_motifs(msa = "msa2.fasta",
                    motifs = ngp)

dplyr::slice(msa1, 120:130)

# A tibble: 11 × 4
   position  seq1  seq2  seq3
      <int> <dbl> <dbl> <dbl>
 1      120     0     0     0
 2      121     0     0     0
 3      122     1     1     1
 4      123     0     0     0
 5      124     0     0     0
 6      125     1     1     0
 7      126     0     0     0
 8      127     0     0     0
 9      128     0     0     1
10      129     0     0     0
11      130     0     0     0

And here's a function to make the plots.
heatmap_separate <- function(data) {
  require(ggplot2)
  require(tidyr)
  numseqs <- ncol(data) - 1
  numrows <- nrow(data)
  toplot <- data |> tidyr::pivot_longer(cols = -position,
                                        names_to = "sequence",
                                        values_to = "hits")
    ggplot(toplot) +
    geom_tile(mapping = aes(x = position,
                            y = sequence,
                            fill = as.factor(hits),
                            color = as.factor(hits))) +
    scale_x_continuous(limits = c(1,numrows),
                       n.breaks = 10,
                       expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0),
                     limits = rev) +
    scale_fill_manual(values = c("#f2f0f7","#54278f")) +
    scale_color_manual(values = c("#f2f0f7","#54278f")) +
    geom_hline(yintercept = 0.5 + 0:numseqs,
               color = "black",
               linewidth = 1.5) +
    theme_minimal() +
    theme(legend.position = "none",
          axis.title.y = element_blank(),
          panel.grid = element_blank(),
          axis.ticks.x = element_line(color = "black"),
          panel.background = element_blank())
}

heatmap_separate(msa1)
heatmap_separate(msa2)

Plot for fasta #1

Plot for fasta #2

The lines representing hits are all plotted as the same color (#54278f), but they appear to be different colors as an artifact of plotting lots of narrow tiles. This can be rectified by viewing the plots as svg files, though I can't display that file type here.


text:
How can I use REMO tool to generate seq.dat and seq.txt files?
I am trying to generate seq.dat and seq.txt files myself, so that I can use them with REMO.
I am using mkseq.pl script in REMOv2.
The following was my command:
user_name@server_name:~/spatial/reconstruction_tools/remo_v2$ perl mkseq.pl /home/spatial/input_pdb/pdb1a6j.pdb 1
PSIPRED finished
user_name@server_name:~/spatial/reconstruction_tools/remo_v2$

This command generated four files:

seq.dat
protein.fasta
psitmp27331007f0101.blast
psitmp27331007f0101.fasta

seq.dat is totally empty. psitmp27331007f0101.blast is also practically empty. There is no seq.txt.
What am I doing incorrectly?

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Error Registering Job Definition AWS
I am in the process of creating a bioinformatics pipeline and I am currently trying to setup the AWS infrastructure for which it will run.
The script I am running to setup the AWS batch job looks like the following:
#! /bin/bash

# Setup the AWS Batch architecture
# Copyright 2022, SolareaBio

set -e

# ***** BEGIN custom lines

container="BGC_bacteria"
vCPUs_req=12
memory_MiB_req=40000
# ASK ABOUT CUSTOM CPU AND MEMORY SETTINGS

# ***** END custom lines

container_lc=$(echo $container | tr [A-Z] [a-z])

subnet_useast1a="subnet-42371d0a"
subnet_useast1b="subnet-b03c54ea"
subnet_useast1c="subnet-00ee1364"

usage() {
  cat << EOF
Usage $0 [-h]

Builds the $container docker container, pushes it to the AWS ECR, and sets up the $container AWS Batch environment
EOF
}

while getopts "h" option; do
  case $option in
    h)   usage;exit;;
    [?]) usage; exit;;
  esac
done

script_dir=$(dirname "$0")

repo="990850896279.dkr.ecr.us-east-1.amazonaws.com"

echo "Logging into ECR"
aws ecr get-login-password --region us-east-1 | docker login --username AWS --password-stdin ${repo}

image_and_tag="${container_lc}:latest"
echo "Building the docker container: $image_and_tag"

# If you get errors from the docker build command, try adding --no-cache
docker build --no-cache -f $script_dir/${container}.dockerfile -t $image_and_tag $script_dir

echo "Tagging Image"
docker tag ${image_and_tag} ${repo}/${image_and_tag}

# Create repository in AWS ECR
repo_container_info=$(aws ecr describe-repositories --repository-names ${container_lc})
if [[ "$repo_container_info" == "" ]]; then
  aws ecr create-repository --repository-name ${container_lc}
fi

echo "Pushing Image up to AWS"
docker push ${repo}/${image_and_tag}

# Create compute environment
EC2_type=ON_DEMAND
ce_info=$(aws batch describe-compute-environments --compute-environments ${container}_${EC2_type}_CE)
ce_status=$(echo $ce_info | jq '.computeEnvironments[0] | .status // ""' | tr -d '"')
if [[ $ce_status != "VALID" ]]; then
  compute_resources='{"type":"EC2","allocationStrategy":"BEST_FIT_PROGRESSIVE","minvCpus":0,"maxvCpus":10000,"tags":{"Name":"'${container}_Batch'"},"instanceTypes":["c5a","m5a","c5","m5"],"subnets":["'$subnet_useast1a'","'$subnet_useast1b'","'$subnet_useast1c'"],"securityGroupIds":["sg-ff2fd48f","sg-65ff5315"],"ec2KeyPair":"SBTools","instanceRole":"ecsInstanceRole"}'
  aws batch create-compute-environment \
    --compute-environment-name ${container}_${EC2_type}_CE \
    --type MANAGED \
    --state ENABLED \
    --compute-resources $compute_resources \
    --service-role AWSBatchServiceRole

  printf "Waiting for Compute Environment to be created"
  while [[ $ce_status != "VALID" ]]; do
sleep 5
printf "."
ce_info=$(aws batch describe-compute-environments --compute-environments ${container}_${EC2_type}_CE)
    ce_status=$(echo $ce_info | jq '.computeEnvironments[0] | .status // ""' | tr -d '"')
  done
  echo ""
fi

# Create the job queue
jq_info=$(aws batch describe-job-queues --job-queues ${container}_${EC2_type}_JQ)
jq_status=$(echo $jq_info | jq '.jobQueues[0] | .status // ""' | tr -d '"')
if [[ $jq_status != "VALID" ]]; then
  aws batch create-job-queue \
    --job-queue-name ${container}_${EC2_type}_JQ \
    --state ENABLED \
    --priority 5 \
    --compute-environment-order '{"order":1,"computeEnvironment":"'${container}_${EC2_type}_CE'"}'
fi

# Register job definition
job_definition="${container}_JD"
jd_info=$(aws batch describe-job-definitions --status ACTIVE --job-definition-name $job_definition)
jd_status=$(echo $jd_info | jq '.jobDefinitions[0] | .status // ""' | tr -d '"')
if [[ $jd_status != "ACTIVE" ]]; then
  aws batch register-job-definition \
    --job-definition-name $job_definition \
    --type container \
    --container-properties '{"privileged":true,"image":"'$repo/$image_and_tag'","vcpus":'${vCPUs_req}',"memory":'${memory_MiB_req}',"mountPoints":[{"sourceVolume":"SB-EFS","containerPath":"/EFS"}],"volumes":[{"host":{"sourcePath":"/"},"name":"SB-EFS"}],"command":["Ref::container","Ref::run_id","Ref::input_file"]}' \
    --retry-strategy '{"attempts": 10,"evaluateOnExit":[{"onStatusReason" :"Host EC2*","action": "RETRY"},{"onReason" : "*","action": "EXIT"}]}'
fi

# Make sure the EFS has an access point set up
efs_id=$(aws efs describe-file-systems | jq -r '.FileSystems[] | select (.Name=="SB-EFS") | .FileSystemId')
if [[ "$efs_id" == "" ]]; then
  echo "ERROR: SB-EFS does not exist"
  exit
fi
ap_id=$(aws efs describe-access-points | jq -r '.AccessPoints[] | select (.Name=="EFS_AP") | .AccessPointId')
if [[ "$ap_id" == "" ]]; then
  # 1001:1001 is user sb in group sb
  ap_id=$(aws efs create-access-point --file-system-id "$efs_id" --posix-user '{"Uid":1001,"Gid":1001}' --root-directory '{"Path":"/"}' --tags '[{"Key":"Name","Value":"EFS_AP"}]' | jq -r '.AccessPointId')
fi

Now when I run it, everything works as it should until getting to the --retry-strategy command line in the "register job definition" section of the script. The following error is subsequently thrown:
Parameter validation failed:
Unknown parameter in retryStrategy: "evaluateOnExit", must be one of: attempts

Not sure what the issue is with the "attempts" parameter that is set. Any potential solutions or guidance to solving this error is much appreciated. Thanks!

Hi @Joe_G, could you please provide a bit more information about what you're trying to do with this script? Stories help a lot with getting good, helpful answers. I can see something suggesting metagenomics, but it's not clear to me what a correct script would do.

My personal approach would be to use the NextFlow AWS pipeline. It would simplify everything, e.g. NextFlow automatically 'Dockerises'

You didn't say what version of the AWS CLI you're using, but I suspect it's version 1. It looks like the evaluateOnExit feature was added in version 2 but was not backported. Installing AWS CLI version 2 resolves the issue:
https://docs.aws.amazon.com/cli/latest/userguide/getting-started-install.html
Note that Nextflow also provides built-in support for AWS Batch. You'd still need to configure some AWS infrastructure using Nextflow, but it should be a lot less painful than writing/maintaining your own script(s). It's definitely worth checking out if you haven't already.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Creating maps.fld file of protein for AutoDock GPU
I want to use AutoDock GPU, which accepts protein.maps.fld and ligand.pdbqt files.
Since I am performing high throughput virtual screening, shall I just load protein AutoDockTools and create protein.maps.fld file, which I can use in my command line run of AutoDock GPU run?
I have 1 million ligands in my screening library, so can't repeat the above procedure on all ligands, or shall I just pick a ligand that contains all types of atoms a ligand can have in that ligand library and produce protein.maps.fld based upon that ligand along with the given protein?


text:
DESeq2 and EdgeR
I am new to using both DESeq2 and EdgeR in Bioconductor used for transforming my RNA expression data.
However, I am struggling to understand their specific purpose, differences between them and underlying theory. The information is required for a thesis with an imminent deadline

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Hi @noor-fatimah, your question seems quite broad and not tied to a specific problem; are you able to make your question more specific, or give some more information on what you've already found out? Have you, for example, read the user's guides for EdgeR and DESeq2?

i need to do proximal vs distal comparison, for that i need to convert my csv file it contains single gene id data multiple time i need to narrow it down to one how can i do that

I think your question has been answered very well here:
https://www.biostars.org/p/284775/
If you want to go deeper perhaps you can look at this paper by Li et al.
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4
where they argue that for large scale datasets (multiple conditions) both of these methods are not optimal and instead suggests using Wilcoxon rank-sum test to decrease chances of FDRs.
Best,
av

i need to do proximal vs distal comparison, for that i need to convert my csv file it contains single gene id data multiple time i need to narrow it down to one how can i do that


text:
How to find closely related genes for a specific gene from WGCNA modules?
I have used WGCNA by integrating several datasets which are processed in a similar way. Identified 17 modules, followed by performed pathway analysis with the genes present in those modules.
Among all the protein-coding genes,ncRNAs and other genes, I'm particularly interested in a specific lncRNA MALAT1.
I found that MALAT1 is co-expressed with more than 800 genes in a module named blue. I have few questions, can anyone please help me in this?
Questions:

I would like to check which are the closest co-expressed genes to MALAT1. How can I know that?

Once I get the closest genes to MALAT1, I would like to make a network with those genes and MALAT1. How to do this network?

The simplest way is to re-calculate (or load, if saved) the TOM matrix and find the genes with the highest TOM to MALAT1. If you have enough RAM to carry out WGCNA in a single block, you can use TOMsimilarityFromExpr, make sure you give it the same arguments as the ones you used in network construction, and select the appropriate column from the result. Then order the column and retain top genes, as many as you think is reasonable. You can also use adjacency instead of TOM, though TOM is a more natural quantity to use if you used in your network construction. Once you figure out which genes you want to retain, you can subset the TOM matrix (both rows and columns) to those genes and you have the sub-network among these genes.

thanks a lot Peter

small question again...I'm interested in modules blue and brown...So, extracted the expression for the genes in those modules and used this way TOM2 = TOMsimilarityFromExpr(t(expr_of_interest),power = softPower,networkType = "unsigned", TOMType = "unsigned") And this gave me correlation values. I can see all correlation values are above 0 only....there are no negative correlation values....why is it like that?

As the function name says, it returns a topological overlap matrix (TOM) which is by definition non-negative. If you want correlation, use cor.

thanks peter. how to get the correlation values for the genes in specific interesting modules or genes from a couple of modules?

Well, you can select the genes using labels, something like (using the notation from WGCNA tutorials) select = moduleColors %in% c(module1, module2,...) (replace module1, module2, .... with the module labels you want to select), then use cor(datExpr[, select]).


text:
Modification time does not trigger a snakemake rule
I'm trying to use snakemake to replace some make workflows, but I'm finding that snakemake does not use the file modification time to trigger a rule.
Here is my snakefile
rule sample:
    input:
        "a"
    output:
        "b"
    shell:
        "touch {output}"

If I create a then run snakemake, b is created as expected:
$ touch a
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job      count    min threads    max threads
-----  -------  -------------  -------------
sample       1              1              1
total        1              1              1

Select jobs to execute...

[Sat Oct 29 14:36:43 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Sat Oct 29 14:36:43 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-29T143643.334050.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 b

Now, if I touch a then ask snakemake to update b, nothing happens:
$ touch a
$ snakemake -c1 b
Building DAG of jobs...
Nothing to be done (all requested files are present and up to date).
Complete log: .snakemake/log/2022-10-29T144005.589943.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 0 Oct 29 14:40 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 b

Modifying the contents of a triggers an update as expected:
$ echo new > a
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job      count    min threads    max threads
-----  -------  -------------  -------------
sample       1              1              1
total        1              1              1

Select jobs to execute...

[Sat Oct 29 14:42:13 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Updated input files: a
    resources: tmpdir=/tmp

[Sat Oct 29 14:42:13 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-29T144213.161985.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 4 Oct 29 14:42 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:42 b

I installed snakemake-minimal (version 7.17.1) from bioconda, am using python 3.9.13 from conda-forge, on a RHEL 7 system.
$ snakemake --version
7.17.1
$ python --version
Python 3.9.13
$ cat /etc/redhat-release
Red Hat Enterprise Linux Server release 7.8 (Maipo)

I'm clearly misunderstanding something very basic about snakemake.  Any suggestions on how to debug this, descriptions of how to have snakemake pay attention to file modification times, or just an explanation of why snakemake behaves differently from make would be much appreciated!
thanks in advance,

This is a Snakemake feature: For small files (less than 100000 bytes), checksums are compared to determine if a job needs to run if the input file is newer than output file:
https://github.com/snakemake/snakemake/pull/1568
With a size of 99999 bytes:
$ dd if=/dev/zero of=a bs=99999 count=1
1+0 records in
1+0 records out
99999 bytes (100 kB, 98 KiB) copied, 0.000502722 s, 199 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:38:00 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Mon Oct 31 14:38:00 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T143800.486966.snakemake.log
$ dd if=/dev/zero of=a bs=99999 count=1
1+0 records in
1+0 records out
99999 bytes (100 kB, 98 KiB) copied, 0.000180857 s, 553 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Nothing to be done (all requested files are present and up to date).
Complete log: .snakemake/log/2022-10-31T143813.285988.snakemake.log

With a size of 100000 bytes:
$ dd if=/dev/zero of=a bs=100000 count=1
1+0 records in
1+0 records out
100000 bytes (100 kB, 98 KiB) copied, 0.000173008 s, 578 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:41:14 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Mon Oct 31 14:41:14 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T144114.211255.snakemake.log
$ dd if=/dev/zero of=a bs=100000 count=1
1+0 records in
1+0 records out
100000 bytes (100 kB, 98 KiB) copied, 5.35361 s, 18.7 kB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:41:26 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Updated input files: a
    resources: tmpdir=/tmp

[Mon Oct 31 14:41:26 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T144125.949828.snakemake.log

Thanks!  Do you know if there is a another mechanism that snakemake uses for "flag files" (https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#flag-files)?

@JulesB Looks like they're just touched or created (see handle_touch()). Then touch() basically just calls os.utime() or touch --no-dereference to set the file access and modified times to the current time.


text:
Running Bionano´s HybridScaffolding from the command line [Error: align0.stdout]
as the title says I´m trying to run Bionanos HybridScaffolding from the command line rather than using their own software interface. This is commercial software to primarily to identify genomewide cancer analysis but also has options for genome assembly.
I installed Perl and all the necessary packages to run the HybridScaffolding Perl script. According to their documentation I started the script with the correct parameters.
After 2 minutes it stops giving me the following error:
Beginning initial NGS CMAP to BioNano CMAP alignment...
Running command: /home/user/workflow/scripts/bionano/tools/pipeline/1.0/RefAligner/1.0/RefAligner
-ref /home/user/results/bionano/fa2cmap/pac_bio_assembly_CTTAAG_0kb_0labels.cmap
-i /home/user/data/bionano/bCalAnn1_Saphyr_DLE1.cmap.gz
-o align0 -stdout -stderr -maxmem 128 -maxthreads 64
-maxvirtmem 0 -maxmemIncrease 8 -RAmem 3 1 0.9 -M 1 3 -ScaleDelta 0.03 2 -ScaleDeltaBPP
-hashScaleDelta 2 -res 2.6 -resSD 0.7 -FP 0.2 -FN 0.02 -sf 0.10 -sd 0.0 -sr 0.01 -se
0.2 -extend 1 -outlier 0.0001 -endoutlier 0.001 -PVendoutlier -deltaX 12 -deltaY 12
-xmapchim 12 -hashgen 5 11 2.4 1.4 0.05 3.0 1 1 1 -hash -hashdelta 27 -hashoffset 1
-hashrange 0 -hashGrouped 6 11 -hashMultiMatch 50 2 -hashGC 300 -hashT2 1 -hashkeys 1
-HSDrange 1.0 -insertThreads 4 -nosplit 2 -biaswt 0 -T 1e-10 -S -1000 -indel -PVres 2
-rres 0.9 -outlierBC -xmapUnique 12 -AlignRes 2. -outlierExtend 12 24 -resEstimate -f
-mres 0.9 -thetaScale 0.5 -KL 1.5 -KF 1.5 -Kmax 4 -MinFN 0.01 -MinFP 0.16 -MinSF 0.05
-MaxSF 0.15 -MinSD 0.0 -MaxSD 0.0 -MinSR 0.007 -MaxSR 0.015 -MaxSE 0.5 -MinRes 2.0
-MinResSD 0.7 -hashSD 0 1 -hashSR 0 1

ERROR: Cannot open align0.stdout: No such file or directory

For context: The HybridScaffolding script contains several scripts which are started automatically over time. One of these is the RefAligner script, which throws the error in this case.
I have a folder called align0 which lies inside the directory that the script is running, therefore I guess this is not the casual "The path is missing" issue..
Inside the Perl Script 'align0.stdout' appears two times:
# error check
if (errCheckRefAligner("align0.stdout", "END of output", "Initial alignment complete in $span.", "ERROR: Initial alignment cannot be completed.", 0) != 1)  {
    my $errMsg = "ERROR: Initial alignment cannot be completed.";
    printStatus(\*STATUS_FILE, "align0", $errMsg, 1);
} else  {
    printStatus(\*STATUS_FILE, "align0", "Done", 0);
} # RefAligner check
errCheckRefAligner("align0.stdout", "END of output", "Initial alignment complete in $span.", "ERROR: Initial alignment cannot be completed.", 1);

Sadly I don´t have any Perl knowledge at all. Does anyone have experience with Bionanos Software and/or knowledge of Perl and could explain to me what is happening inside this part of the script?

What operating system are you using? Are you sure the script can run on it? The error suggests perl is trying to open a file named align0.stdout for either reading or writing and that file doesn't exist. Since you are passing align0 as an output file name with -o align0 (I am guessing), then presumably the script will create $outputName.stdout somewhere. Does it work if you give a full path like -o /home/user/data/bionano/align0? Also, what directory are you running this from? D you have write permissions there?

I´m working on a Debian Server. I was able to spot the problem. I thought the Server had at least 256 RAM. Turns out it had only 64. I think HybridScaffolding needs at least 128 RAM. So when I ran the whole HybridScaffolding Script on a stronger Server (1000 RAM) it worked just fine without any problems :)

I was able to spot the problem. I thought the Server I was working on had at least 256 RAM. Turns out it only had 64.
I think HybridScaffolding needs at least 128 RAM. So when I ran the whole HybridScaffolding script on a stronger Server (1000 RAM) it worked just fine without any problems :)
Edit: For the record, I also had to unzip the cmap.gz file because the RefAligner could not read it zipped.


text:
How to count number of motifs in Prosite?
For example I checked motifs in Transferrin P02787  https://prosite.expasy.org/cgi-bin/prosite/ScanView.cgi?scanfile=816050634835.scan.gz
I have
hits by profiles: [2 hits (by 1 profile) on 1 sequence]
and
hits by patterns: [5 hits (by 3 distinct patterns) on 1 sequence]
Can I say that I detect 4 different motifs? One motif is searched by profile searching (PS51408   TRANSFERRIN_LIKE_4   Transferrin-like domain profile two times) and three motifs are searched by pattern searching (PS00205   TRANSFERRIN_LIKE_1   Transferrin-like domain signature 1, PS00206   TRANSFERRIN_LIKE_2   Transferrin-like domain signature 2 two times, PS00207   TRANSFERRIN_LIKE_3   Transferrin-like domain signature 3 two times).
All the motifs found are equal to 7? Is it right? Can I talk about the different motifs found by the profile and pattern, or are they the same in some way?

As per PROSITE database manual:

Unlike patterns, profiles are usually not confined to small regions
with high sequence similarity. Rather they attempt to characterize a
protein family or domain over its entire length. This can lead to
specific problems not arising with PROSITE patterns. With a profile
covering conserved as well as divergent sequence regions, there is a
chance to obtain a significant similarity score even with a partially
incorrect alignment. This possibility is taken into account by our
quality evaluation procedures. In order to be acceptable, a profile
must not only assign high similarity scores to true motif occurrences
and low scores to false matches. In addition, it should correctly
align those residues having analogous functions or structural
properties according to experimental data.
Profiles are supposed to be more sensitive and more robust than
patterns because they provide discriminatory weights not only for the
residues already found at a given position of a motif but also for
those not yet found. The weights for those not yet found are
extrapolated from the observed amino acid compositions using empiric
knowledge about amino acid substitutability.

Given the less biased and somewhat more balanced result provided by a profile I would say that the sequence contains two motifs that correspond to a PS51408 TRANSFERRIN_LIKE_4 Transferrin-like domain profile and provide the bitscore and e-value for the match.
Best,
av


text:
How to adjust #columns shown when PyRange is printed?
In Pandas, we have pandas.options.display to adjust how a DataFrame should be displayed when printed.
I wonder if PyRanges has a similar set of options?
For example, if I have a PyRange object like this:
import pyranges as pr
data = pr.data.ucsc_bed()

I can do print(data) and get:
+--------------+-----------+-----------+------------+------------+-------+
| Chromosome   | Start     | End       | Feature    | gene_id    | +4    |
| (category)   | (int32)   | (int32)   | (object)   | (object)   | ...   |
|--------------+-----------+-----------+------------+------------+-------|
| chr1         | 12776117  | 12788726  | gene       | AADACL3    | ...   |
| chr1         | 169075927 | 169101957 | gene       | ATP1B1     | ...   |
| chr1         | 6845383   | 7829766   | gene       | CAMTA1     | ...   |
| chr1         | 20915589  | 20945396  | gene       | CDA        | ...   |
| ...          | ...       | ...       | ...        | ...        | ...   |
| chrX         | 152661096 | 152663330 | exon       | PNMA6E     | ...   |
| chrX         | 152661096 | 152666808 | transcript | PNMA6E     | ...   |
| chrX         | 152664164 | 152664378 | exon       | PNMA6E     | ...   |
| chrX         | 152666701 | 152666808 | exon       | PNMA6E     | ...   |
+--------------+-----------+-----------+------------+------------+-------+
Stranded PyRanges object has 5,519 rows and 9 columns from 30 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.
4 hidden columns: transcript_id, Strand, exon_number, transcript_name

How can I increase the number of columns shown?
Also appreciated if you refer me to a documentation that explains other related display options.
Update: For future reference, @terdon remark encouraged me to check if such an adjustment is even considered in the PyRanges. I found the line that is responsible for this representation, and indeed the COLUMN environment variable is the only defining parameter. In conclusion, the adjustment that I am looking for is simply not implemented in PyRanges.

Are you running in an interactive python shell? On what operating system? What happens if you make the terminal window larger? Or if you print to a file? Can you provide a minimal code example that we can run which reproduces the problem?

Thanks @terdon for your reply. I already provided the minimal code example. Just run the code I provided and you get this picture.

The only code you have provided is data = pr.data.ucsc_bed(). We would need an example file, and the commands to read it into the script etc. Please see https://stackoverflow.com/help/minimal-reproducible-example. Also, please answer my other questions too.

@terdon, please note that ucsc_bed() produces a dummy data. I now added import pyranges as pr. Just copy-paste this code and you get the output. And regarding your questions: I am answering :). Sorry, I pressed Enter and the comment got posted. I tested this in the terminal and also jupyter notebook, but the output is independent of the environment. It also stays the same irrespective of the size of window.

Ah! Thank you, I hadn't realized that. And thanks for the answers! Does it print the whole thing if printed to a file, at least?

Of course. You can always get the Pandas.DataFrame using data.df, and then display it as you wish. Alternatively, you can export the bed file using data.to_bed('./test.bed'). In all cases, all data is exported. Note that my question is related to display of the data in the terminal. Not exporting issues. Please consider taking a look at pandas.options.display to understand what I mean by difference between exporting the PyRanges objects, and displaying PyRanges objects.

Let us continue this discussion in chat.

When trying this on my machine, if I resize the terminal window before running print(data), I get the full output:

You can also force this, even if your terminal is actually smaller, by setting the COLUMNS variable to something larger than the actual size when layching python or ipython. This is a special bash variable that holds the current size of the terminal. From man bash:
   COLUMNS
          Used  by  the select compound command to determine the terminal
          width when printing selection lists.  Automatically set if  the
          checkwinsize  option is enabled or in an interactive shell upon
          receipt of a SIGWINCH.

So if you set it manually, you can get the entire data frame printed, BUT you don't really want to do that since it won't scroll, it just wraps the text. For example, try:
COLUMNS=300 python

That will set the size to 300, which results in this:

Very good @terdon. Thanks! I tested this with Jupyter notebook and even when adjusting the window size following this the additional columns are still kept hidden. Most likely, I should change another parameter in Jupyter. But indeed, I could confirm that in Terminal if I adjust the window size (even inside a screen session), the columns do show as long as the window size is large enough. Your approach using COLUMN variable is very nice!

You're welcome. If Jupyter gives you a shell, you can try setting the variable in the shell. Alternatively, you can actually do it in your python code with os.environ['COLUMNS']="300" and then run print(data) and it might work as expected (confirmed in an ipython session, but I don't have a jupyter notebook to test on).

I confirm that setting the environment variable using os.environ['COLUMNS']="300" works in a Jupyter notebook too! Thank you very much @terdon.


text:
How to identify matching and unique rows across different columns in R data-frame?
I am working on a large data-frame that has 4 columns and each column has variable rows. I am trying to find unique and identical pathways between the 4 columns (each column represents a particular day of treatment with a drug). Here below is a small example.
Pathways_1Day <- c("blood","kidney","testis","No","bone","liver","intestine","lungs","ABC","pancreas","Yes","spleen")

Pathways_2Day <- c("blood","kidney","testis","eyes","bone","cells","intestine","cervix","ABC","pancreas","None")

Pathways_3Day <- c("blood","kidney","vessels","lymph","t-cells","liver","intestine","lungs","ABC","epidermis","None")

df<-data.frame(Pathways_1Day,Pathways_2Day,Pathways_3Day)

I want to get a summary of the no of pathways that are common between different timepoints (1, 2 and 3 days).
Important: The no of pathways is not the same for each day.
I have tried this:
All_pathwayNames <- df%>%group_by_all%>%count

But the desired output is not what I am trying to get.
There can be different ways to address that. It will be great if I can get matching rows infront of each other across all columns.
Regards

Are pathways on different day with the same genes? Then It is just a matter to compare the names of each day via setdiff(Pathways_3Day, Pathways_1Day) for example. If not I would recommend it to work with lists instead of data.frame for this kind of data or use GeneSetCollection from GSEABase to read it and work with this.

When I create the dataframe using the commands you have provided, I get an error, and the data frame can't be created:
Error in data.frame(Pathways_1Day, Pathways_2Day, Pathways_3Day) :
  arguments imply differing number of rows: 12, 11

What I think you were trying to do is something more like this:
library(tidyverse)
tibble() %>%
  bind_rows(tibble(day=1, pathway = Pathways_1Day)) %>%
  bind_rows(tibble(day=2, pathway = Pathways_2Day)) %>%
  bind_rows(tibble(day=3, pathway = Pathways_3Day)) -> combined.tbl

Producing a long-format data structure that looks like this:
> combined.tbl
# A tibble: 34 x 2
     day pathway
   <dbl> <chr>
 1     1 blood
 2     1 kidney
 3     1 testis
 4     1 No
 5     1 bone
 6     1 liver
 7     1 intestine
 8     1 lungs
 9     1 ABC
10     1 pancreas
# ... with 24 more rows

When laid out in this way, you can do something like grouping by pathway, and then do a pivot_wider to show the sets:
> combined.tbl %>%
     group_by(pathway) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 20 x 4
# Groups:   pathway [20]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 No            1    NA    NA
 3 None         NA     2     3
 4 Yes           1    NA    NA
 5 blood         1     2     3
 6 bone          1     2    NA
 7 cells        NA     2    NA
 8 cervix       NA     2    NA
 9 epidermis    NA    NA     3
10 eyes         NA     2    NA
11 intestine     1     2     3
12 kidney        1     2     3
13 liver         1    NA     3
14 lungs         1    NA     3
15 lymph        NA    NA     3
16 pancreas      1     2    NA
17 spleen        1    NA    NA
18 t-cells      NA    NA     3
19 testis        1     2    NA
20 vessels      NA    NA     3

And insert a filter to select only pathways that appear in all three days (i.e. n() == 3):
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() == 3) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 4 x 4
# Groups:   pathway [4]
  pathway     `1`   `2`   `3`
  <chr>     <dbl> <dbl> <dbl>
1 ABC           1     2     3
2 blood         1     2     3
3 intestine     1     2     3
4 kidney        1     2     3

Or at least two days (n() >= 2):
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() >= 2) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 10 x 4
# Groups:   pathway [10]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 None         NA     2     3
 3 blood         1     2     3
 4 bone          1     2    NA
 5 intestine     1     2     3
 6 kidney        1     2     3
 7 liver         1    NA     3
 8 lungs         1    NA     3
 9 pancreas      1     2    NA
10 testis        1     2    NA

You're misusing data frames. Split the df into 4 vectors and use a Venn diagram generator to denerate a Venn diagram with the 4 vectors.
Adapting the first example from here
#untested
venn.diagram(
  x = list(df[,1], df[,2], df[,3], df[,4]),
  category.names = colnames(df),
  filename = 'my_venn_diagram.png',
  output=TRUE
)

Good call @RamRS


text:
Why does gene abundance table from StringTie output the same gene ID twice?
I would like to compare gene level counts (FPKM) from StringTie output. I understand FPKM is outdated but my PI prefers to use it as a reference/guide in conjunction with the normalised counts from DESeq2.
I processed my samples by first running StringTie on my aligned read files (*.bam) with reference to the ensemble GTF:
stringtie sample$_{i}.sorted.bam -p 12 -G reference.gtf --rf -o $OUT_DIR/sample_${i}/sample_${i}.gtf -v -B

I left out the -e flag in order for novel transcripts to be assembled.
I then merged the individual GTFs using:
stringtie --merge -G referece.gtf -o merged.gtf -p 12 -i assembly_GTF_list.txt

Finally I re-run StringTie using the merged.gtf as reference to obtain expression values:
stringtie sample_${i}.sorted.bam -p 12 -G merged.gtf --rf -o sample_{i}.gtf -v -e -B -A sample_{i}.gene_abund.tab

However, in the gene_abund.tab of one sample, I find that the same Gene ID is quantified twice (albeit with different coordinates):
Gene ID Gene Name       Reference       Strand  Start   End     Coverage        FPKM    TPM
ENSG00000225657 LINC01716       20      +       56577563        56596126        0.000000        0.000000        0.000000
ENSG00000225657 LINC01716       20      +       56596310        56598621        0.0     0.0     0.0

These are the transcripts corresponding to 'ENSG00000225657' in the reference gtf I used:
20      havana  transcript      56577563        56596126        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; gene_name "LINC01716"; ref_gene_id "ENSG00000225657";
20      havana  exon    56577563        56577612        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "1"; gene_name "LINC01716";
20      havana  exon    56577742        56577950        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "2"; gene_name "LINC01716";
20      havana  exon    56579391        56579553        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "3"; gene_name "LINC01716";
20      havana  exon    56580201        56580335        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "4"; gene_name "LINC01716";
20      havana  exon    56595063        56595130        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "5"; gene_name "LINC01716";
20      havana  exon    56596054        56596126        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000434176"; exon_number "6"; gene_name "LINC01716";
20      havana  transcript      56596310        56598621        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000437350"; gene_name "LINC01716"; ref_gene_id "ENSG00000225657";
20      havana  exon    56596310        56596488        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000437350"; exon_number "1"; gene_name "LINC01716";
20      havana  exon    56598503        56598621        .       +       .       gene_id "ENSG00000225657"; transcript_id "ENST00000437350"; exon_number "2"; gene_name "LINC01716";

I assumed that gene level counts would be the aggregation of all transcript level counts. However, it looks like StringTie is reporting the transcript level counts instead. Have I missed out a step here? Or is the gene_abund.tab not meant to contain gene-level counts?

You have two transcripts for that gene, so presumably the output is one line per transcript, right? If you really want to treat them as one things, can't you just sum up the values?

@terdon yes, I guess I could sum them. However, I was expecting these to be gene level counts so I'm not sure if I've messed up something in the workflow.


text:
How to find restriction enzyme frequency in whole genome and count the distance between them?
I am having possibly a simple question which is proving difficult.
My objective is to take the reference FASTA file and flat file database with 1500 restriction enzymes with its cutting sites (mostly pattern I would say).

I first want to see that how many times each enzyme cut the FASTA file and;
second to check the number and size of fragment that each enzyme would give me, this output pregerably would be a CSV file may be.

Background I tried to do pattern search using R, but its just collapsing my system. I am looking for bash command may be.

This looks like the RAM load of the reg-ex in R. Bash isn't the best for reg-ex style commands IMO, Python can do this, I would assume awk too, and dare I say Perl is really good.

Biopython has a module dedicated to working with restriction enzymes. Without any code examples we cannot help you much further.

The diffHic package might have solved this issue in R.
# Load the genomic data
genome <- biostrings::readDNAStringSet("your_file.fa")

# Run the in silico restriction digest
cuts <- cutGenome(genome, "CCGCGG", overhang=2L) # SacII

the cuts GRanges object contains the boundaries of each restriction fragment for given genome.


text:
Using SED command in Nextflow script
I was wondering if someone can help me figure out this error message. I’m at the end of my Nextflow pipeline and I want to change the header in a FASTA file and it works when I use this command:
sed -i "s/^>/>${assembly_id}_${gene_id}_/g" "${assembly_id}/${gene_id}/${assembly_id}_${gene_id}.fasta"
but it prints out much more than I actually need (i.e., “goodID_extra_extra_extra”). That’s just part of the file name (it includes what assembler was used) but I only want it to add “goodID_” to the header.
I tried using this command:
sed -i "s/^>/>${assembly_id%%_*}_${gene_id}_/g" "${assembly_id}/${gene_id}/${assembly_id}_${gene_id}.fasta"
but I get this error message:
N E X T F L O W  ~  version 21.10.6
Launching `geneTargets.nf` [amazing_elion] - revision: b875027a52
Script compilation error
- file : /Users/rimo/Desktop/AWS/nextflow_work/PROJECTS/gene_targets/geneTargets.nf
- cause: Unexpected input: '{' @ line 94, column 32.
   process publish_gene_sequences {
                                  ^

1 error

Can someone explain why this is happening? If this is something I can’t do without breaking the script that’s okay the header will still do it’s purpose I just wanted it to be clearer/less noise.
EDIT: this is the process block
process publish_gene_sequences {
publishDir params.publish_dir
input:
    tuple val(assembly_id), val(gene_id), path("gene.fasta")
output:
    path("${gene_id}/${assembly_id}_${gene_id}.fasta")
script:
    """
    mkdir -p ${gene_id}
    cp gene.fasta "${gene_id}/${assembly_id}_${gene_id}.fasta"
    sed -i "s/^>/>${assembly_id}_${gene_id}_/g" "${gene_id}/${assembly_id}_${gene_id}.fasta"
    """
 }

Can you please share a minimal reproducible example? For now, I would guess this has something to do with escaping like this example here: https://github.com/nextflow-io/nextflow/issues/2082#issuecomment-1025569003

I added the entire process block where the line of code I'm trying to use is @mribeirodantas

Ideally, a minimal reproducible example would include fold structure, some sample files and a workflow call to this process block. Trying to guess things on my side, adding -e to the sed line after -i gave me no errors but it's guess work unless I can try to reproduce what you're finding in your machine :-\

Inside your script's command string, assembly_id is actually a java.lang.String so trying to manipulate it with a Bash operation (i.e. using ${assembly_id%%_*}) won't produce the expected result. Here's one way to return a substring up to the first underscore; note that we can use single quotes around our sed expression:
process publish_gene_sequences {

    publishDir params.publish_dir

    input:
    tuple val(assembly_id), val(gene_id), path("gene.fasta")

    output:
    path("${gene_id}/${assembly_id}_${gene_id}.fasta")

    script:
    def assembly = assembly_id.substring(0, assembly_id.indexOf('_'))

    """
    mkdir "${gene_id}"
    cp gene.fasta "${gene_id}/${assembly_id}_${gene_id}.fasta"
    sed -i 's/^>/>${assembly}_${gene_id}_/g' "${gene_id}/${assembly_id}_${gene_id}.fasta"
    """
}

Note also that in the above example, copying the file can be avoided by applying the expression directly and redirecting stdout:
script:
def assembly = assembly_id.substring(0, assembly_id.indexOf('_'))

"""
mkdir "${gene_id}"
sed \\
    -e 's/^>/>${assembly}_${gene_id}_/g' \\
    gene.fasta \\
    > "${gene_id}/${assembly_id}_${gene_id}.fasta"
"""

That makes a lot more sense thank you for clarifying! That solution worked for me!


text:
Creating a PDB from a figure in paper
I want to create a PDB file of a drug that is given in a paper.
Can somebody share a tutorial or similar resources for a beginner-level user?

Thank you.

Could you add some context? What do you want the PDB file for? Do you have the structure of this molecule in another file format? (I didn't read the paper).

yep have same question. Wouldn' t a cif file better serve the use if needed for fitting a ligand into electron density ? Or you just nred to depict molecule in 3D ?

You can use Image2Smiles tool. It produces smi format, then you can convert smi to pdb using openbabel.
If you just have a few figures, I suggest using PubChem Sketcher.

or you could taket it from :
https://www.rcsb.org/structure/7WMV pdb file
HETATM 4865  C13 1YI A 701     106.585 114.218 112.698  1.00 74.13           C
HETATM 4866  C15 1YI A 701     106.050 112.633 114.547  1.00102.13           C
HETATM 4867  C17 1YI A 701     104.024 111.140 110.766  1.00 89.07           C
HETATM 4868  C20 1YI A 701     103.626 108.779 110.527  1.00 83.00           C
HETATM 4869  C21 1YI A 701     104.626 108.507 111.437  1.00 84.86           C
HETATM 4870  C22 1YI A 701     102.234 110.466 109.205  1.00 85.32           C
HETATM 4871  C24 1YI A 701     100.911 110.073 109.846  1.00 66.83           C
HETATM 4872  C26 1YI A 701      92.545 108.407 114.139  1.00 95.68           C
HETATM 4873  C28 1YI A 701     100.141 109.053 109.319  1.00 63.43           C
HETATM 4874  C11 1YI A 701     105.694 113.343 111.806  1.00 75.91           C
HETATM 4875  C12 1YI A 701     105.840 111.929 112.354  1.00 91.76           C
HETATM 4876  C14 1YI A 701     105.990 114.087 114.099  1.00 85.43           C
HETATM 4877  C16 1YI A 701     105.029 110.853 111.663  1.00 99.03           C
HETATM 4878  C18 1YI A 701     105.319 109.548 112.006  1.00100.47           C
HETATM 4879  C19 1YI A 701     103.333 110.096 110.188  1.00 81.96           C
HETATM 4880  C23 1YI A 701     104.801 110.649 116.200  1.00 88.38           C
HETATM 4881  C25 1YI A 701     102.892 107.590 109.921  1.00 99.79           C
HETATM 4882  C27 1YI A 701     100.517 110.745 110.979  1.00 88.80           C
HETATM 4883  C29 1YI A 701      98.556 109.383 111.080  1.00 77.66           C
HETATM 4884  C30 1YI A 701      97.245 109.006 111.756  1.00 83.68           C
HETATM 4885  C31 1YI A 701      97.108 107.503 111.935  1.00 85.27           C
HETATM 4886  C32 1YI A 701      99.328 110.410 111.602  1.00 90.75           C
HETATM 4887  C33 1YI A 701      98.954 108.713 109.937  1.00 76.18           C
HETATM 4888  C34 1YI A 701      96.107 107.127 113.039  1.00 98.03           C
HETATM 4889  C35 1YI A 701      94.780 107.913 113.042  1.00 92.89           C
HETATM 4890  C36 1YI A 701      91.696 107.883 115.298  1.00113.05           C
HETATM 4891  C37 1YI A 701      92.889 109.868 114.425  1.00 87.90           C
HETATM 4892  C38 1YI A 701      91.711 108.277 112.861  1.00 77.85           C
HETATM 4893  N8  1YI A 701      93.804 107.673 114.094  1.00 93.06           N
HETATM 4894  N9  1YI A 701      90.968 106.633 115.149  1.00106.00           N
HETATM 4895  O2  1YI A 701     105.306 111.885 113.648  1.00 99.43           O
HETATM 4896  O3  1YI A 701     106.050 113.444 110.457  1.00 80.66           O
HETATM 4897  O4  1YI A 701     106.560 115.552 112.284  1.00 86.11           O
HETATM 4898  O5  1YI A 701     106.634 114.878 115.050  1.00 99.48           O
HETATM 4899  O6  1YI A 701      94.512 108.703 112.208  1.00 84.41           O
HETATM 4900  O7  1YI A 701      91.625 108.487 116.316  1.00101.67           O
HETATM 4901  S1  1YI A 701     105.267 112.404 116.160  1.00121.62           S

picture :

obtained with : https://litemol.org/Viewer/
OOOPSSSSSSSSSS
realized that ligand in pdb its not the full ligand of your picture,
need to get it from Ligand page on RCSB_PDB : https://www.rcsb.org/ligand/1YI
download links highlighted in pics :

download /1YI_ideal.sdf  convert via openbabel :
http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html
result
WARNING: damped steplength 0.0184262 to 0.00725371
WARNING: damped steplength 0.00508674 to 0.00100138
COMPND    1YI
AUTHOR    GENERATED BY OPEN BABEL 2.3.2
HETATM    1  C   LIG     1       3.537  -2.101   2.235  1.00  0.00           C
HETATM    2  C   LIG     1       2.812  -4.346   1.266  1.00  0.00           C
HETATM    3  C   LIG     1       2.262  -2.902  -2.712  1.00  0.00           C
HETATM    4  C   LIG     1       0.447  -1.602  -3.663  1.00  0.00           C
HETATM    5  C   LIG     1       0.405  -0.869  -2.474  1.00  0.00           C
HETATM    6  C   LIG     1       1.522  -3.533  -5.065  1.00  0.00           C
HETATM    7  C   LIG     1       0.305  -3.796  -5.950  1.00  0.00           C
HETATM    8  C   LIG     1      -5.920  -1.833 -12.368  1.00  0.00           C
HETATM    9  C   LIG     1      -0.581  -4.848  -5.684  1.00  0.00           C
HETATM   10  C   LIG     1       2.908  -1.583   0.938  1.00  0.00           C
HETATM   11  C   LIG     1       3.139  -2.484  -0.313  1.00  0.00           C
HETATM   12  C   LIG     1       3.714  -3.615   2.273  1.00  0.00           C
HETATM   13  C   LIG     1       2.202  -2.172  -1.496  1.00  0.00           C
HETATM   14  C   LIG     1       1.247  -1.156  -1.414  1.00  0.00           C
HETATM   15  C   LIG     1       1.366  -2.666  -3.802  1.00  0.00           C
HETATM   16  C   LIG     1       1.091  -5.111   3.310  1.00  0.00           C
HETATM   17  C   LIG     1      -0.531  -1.162  -4.707  1.00  0.00           C
HETATM   18  C   LIG     1      -0.009  -2.939  -7.020  1.00  0.00           C
HETATM   19  C   LIG     1      -2.131  -4.074  -7.418  1.00  0.00           C
HETATM   20  C   LIG     1      -3.479  -4.118  -8.115  1.00  0.00           C
HETATM   21  C   LIG     1      -3.492  -4.837  -9.473  1.00  0.00           C
HETATM   22  C   LIG     1      -1.224  -3.051  -7.715  1.00  0.00           C
HETATM   23  C   LIG     1      -1.778  -4.991  -6.415  1.00  0.00           C
HETATM   24  C   LIG     1      -4.883  -4.859 -10.129  1.00  0.00           C
HETATM   25  C   LIG     1      -5.503  -3.491 -10.484  1.00  0.00           C
HETATM   26  C   LIG     1      -7.501  -1.753 -12.406  1.00  0.00           C
HETATM   27  C   LIG     1      -5.427  -1.719 -13.834  1.00  0.00           C
HETATM   28  C   LIG     1      -5.340  -0.600 -11.644  1.00  0.00           C
HETATM   29  N   LIG     1      -5.408  -3.127 -11.821  1.00  0.00           N
HETATM   30  N   LIG     1      -8.129  -2.821 -13.038  1.00  0.00           N
HETATM   31  O   LIG     1       3.072  -3.897  -0.073  1.00  0.00           O
HETATM   32  O   LIG     1       3.421  -0.244   0.755  1.00  0.00           O
HETATM   33  O   LIG     1       2.698  -1.683   3.334  1.00  0.00           O
HETATM   34  O   LIG     1       5.094  -3.922   1.972  1.00  0.00           O
HETATM   35  O   LIG     1      -6.061  -2.821  -9.650  1.00  0.00           O
HETATM   36  O   LIG     1      -8.120  -0.753 -12.076  1.00  0.00           O
HETATM   37  S   LIG     1       1.026  -4.388   1.646  1.00  0.00           S
HETATM   38  C   LIG     1      -9.774  -4.156 -16.922  1.00  0.00           C
HETATM   39  C   LIG     1      -9.839  -1.879 -16.189  1.00  0.00           C
HETATM   40  N   LIG     1      -9.383  -3.231 -15.850  1.00  0.00           N
HETATM   41  C   LIG     1      -9.964  -3.665 -14.556  1.00  0.00           C
HETATM   42  C   LIG     1      -9.568  -2.827 -13.310  1.00  0.00           C
CONECT    1   10   12   33
CONECT    2   31   12   37
CONECT    3   15   13
CONECT    4   17   15    5
CONECT    5    4   14
CONECT    6    7   15
CONECT    7    6    9   18
CONECT    8   28   29   27   26
CONECT    8
CONECT    9    7   23
CONECT   10   11    1   32
CONECT   11   13   10   31
CONECT   12    1    2   34
CONECT   13    3   14   11
CONECT   14    5   13
CONECT   15    6    4    3
CONECT   16   37
CONECT   17    4
CONECT   18    7   22
CONECT   19   23   22   20
CONECT   20   19   21
CONECT   21   20   24
CONECT   22   18   19
CONECT   23    9   19
CONECT   24   21   25
CONECT   25   35   24   29
CONECT   26    8   30   36
CONECT   27    8
CONECT   28    8
CONECT   29   25    8
CONECT   30   26   42
CONECT   31   11    2
CONECT   32   10
CONECT   33    1
CONECT   34   12
CONECT   35   25
CONECT   36   26
CONECT   37    2   16
CONECT   38   40
CONECT   39   40
CONECT   40   38   39   41
CONECT   41   40   42
CONECT   42   41   30
MASTER        0    0    0    0    0    0    0    0   42    0   42    0
END

picture here :

more about the compound at:
PDB Chemical Component 1YI
UPDATE
Last link doesnt work anymore :

Ligand Expo
Ligand Expo is down from January 1-31, 2023. Requests to this service will return a 410 http error code (“Gone”). Please help us better understand which Ligand Expo features are utilized and should be integrated with RCSB.org by taking this survey.

Users can access CCD data in mmCIF format at https://files.wwpdb.org/pub/pdb/refdata/chem_comp/. For other data download during server downtime, please contact chenghua.shao@rcsb.org.

Can use : https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/1YI


text:
In CellRanger output, what are the reads with an xf:i:17 tag?
When running CellRanger on 10x scRNA-Seq data, a BAM file is output with the aligned reads. These reads have an xf tag, described here:

The bits of this tag are interpreted as follows:

1 The read is confidently mapped to the transcriptome
2 This read's barcode, UMI, and feature combination was discarded in favor of a different feature with higher read support
4 This read pair maps to a discordant pair of genes, and is not treated as a UMI count
8 This read is representative for a molecule and is treated as a UMI count
16 This read maps to exactly one feature, and is identical to bit 1 for transcriptomic reads. Notably, this bit is set for a Feature
Barcode read, while bit 1 is not
32 This read was removed by targeted UMI filtering.

In a real dataset, I have this distribution of xf values:
$ samtools view possorted_genome_bam.bam |  egrep -o "xf:i:[0-9]+" | sort | uniq -c
66263599 xf:i:0
99892363 xf:i:17
  245179 xf:i:19
 8283967 xf:i:25

It is clear to me that 0 are unaligned, 19 = 16+2+1 are discarded, 25 = 16+8+1 are kept (and used in the feature_barcode_matrix); but I do not understand what happens to the 17 = 16 +1. They seem to be good reads, they are aligned on genes, yet they are discarded.
One possible explanation from my reading of this paragraph is that the reads with 17 are PCR duplicates, thus ignored, whereas the reads with 19 are in a group matching different genes (several reads with the same BC/UMI but matching different genes). But this doesn't look right when exploring the BAM file manually.

From the documentation at https://kb.10xgenomics.com/hc/en-us/articles/115003646912-How-is-sequencing-saturation-calculated-

All of these samflag 0x400 reads have an xf tag value of 17, which consist of bits 1 and 16. This also means these reads do not have the xf bit of 8, which mark representative reads from a group of duplicates. The converse isn't true though. The xf17 consist mostly of samflag 0x400 duplicate reads but also of samflag nonduplicate reads.

Thank you for the useful answer, however I don't think it's fully answering the question: that suggests that some of these reads, but not all, are duplicates of some sort, but still leaves to me some mystery about what the other ones are, how the duplicates are detected etc


text:
What are the columns in bedtools coverage output hist file?
I am using bedtools to caculate the coverage of my targets of my WES data, to later plot. But to plot, I need to know what to plot and what it is I am seeing. I have unsuccefully tried to find what exactly the columns are in the output of the bedtool coverage, and how to read them. Also, the documentation does not help me a lot, and I find it confusing. Does anyone know, or where to find this a bit more elaboratly explained?
I used the following code:
bedtools coverage -hist -abam input.bam -b input.bed > output >2 log
Correction note: the code above is actually wrong (thanks to @Alexlok) and should be:
bedtools coverage -hist -a input.bed -b input.bam > output 2> log
And the output is (first line):

Last 5 lines:

Thank you in advance

Make sure that you're not following a script written for bedtools version < 2.24: the inputs a and b are inverted in recent versions. So on a recent version of bedtools, if you want to count the read coverage of your features, you need: bedtools coverage -a input.bed -b input.bam. Your current command counts the number of features in b that intersect each read in a, and I'm not sure that what you want to do.

I am following workflows for exome coverage probally indeed using <2.24, since they were published a while ago... not finding any new versions. I want to see how much of the exome targets (bed file) are actually captured by my reads indeed. I use -abam (my aligned reads to the reference genome) and -b (my target bed file) option? So I am not actually using the -a option. Does it also count then? Thanks

If I understand correctly (I have only used bedtools a bit, no expert on that), on older versions (before 2.19.0) -a would have expected a bed file, so -abam made explicit that it was a bam. But on more recent versions, -a and -b automatically detect the input type, so -abam is synonym with -a (kept for compatibility). And since version 2.24 the options a and b are inverted, so you would need a -bbam (which does not exist since -b covers that use case). See the example in the page you linked.

Thank you so much for that. You are absolutly right and I got totally confused.

Found it, it is an extended .bed file format, so not specific to bedtools output: source bed file output.
With 4 added columns to the end from bedtools -hist option:
bedtools -hist option

Thanks for posting the answer, appreciated.


text:
Can I discard missing alleles?
I am converting a biallelic VCF into a SQL table and one of my tables will be something like:
SAMPLE | custom variant identifier | Genotype
-------|---------------------------|---------
0001   | 3456                      | ./.

My question is that for the diploid ./. or haploid . missing alleles, whether I can just drop them to save space. My thinking is that since the dot represents missing data, the ./. are just placeholders in the VCF but provide no information at all.
Is this correct? Should I retain the rows with ./.?
Showing the end of one line of the VCF (really long lines because of high N):
FORMAT          SampleX             SampleY
GT:AD:DP:GQ:PL  0/0:3,0:3:4:0,4,45  ./.:3,0:3:.:0,0,0

I thought to also show the result after bcftools norm -m to split to biallelic and queried per sample.
sample_id   chr@pos@ref@alt genotype
sample1 chr1@10158@ACCCT@A       0/0
sample2 chr1@10158@ACCCT@A       ./. <- does this provide any info or can I discard?
sample3 chr1@10158@ACCCT@A       0/0
sample4 chr1@10158@ACCCT@A       ./.
sample5 chr1@10158@ACCCT@A       0/0

Can you show us a few lines of the actual VCF file please? ./. means the site wasn't genotyped. It is common in gvcf files, for example but we'd have to see some data to judge correctly. Where are the "asterisks"? Do you mean .? Asterisks are *.

Thanks, I've updated the question. I guess what I'm asking is whether there is any information being encoded for the samples that weren't genotyped for that variant i.e. ./.

If you can, do not discard missing allele information.
Generally you often benefit from being able to distinguish between a position that you interrogated and genotyped as REF (reference) versus one that you interrogated but could not determine the genotype (e.g due to low coverage in the sequencing or probing issues in a microarray setup) and which usually is encoded as a missing (./. if ploidy = 2) genotype. A missing genotype tell us that the genotype at this position could be REF, but could also be some other base, we simply do not know.
Knowing this is critical and/or useful depending the downstream analysis. For example it is quite important and informative when merging datasets. Since formats like VCF or PLINK do not record non-variant positions, if you discarded missing genotypes prior to merging, how can you tell between genotypes that truly were REF in one dataset (fully invariant position) and some recorded missing genotypes at the same position but which you decided to discard? You will end up assuming REF status for everything not recorded and this will inevitably lead to batch effects in your combined dataset. For a better discussion see this blog post
On the other hand, if your relational database (your particular situation) for instance will not incorporate external databases in order to produce joint datasets, perhaps discarding missing genotypes is all fine.


text:
Removing particular line along with Fasta header in fasta file
I have a fasta file of 3'UTR sequences of the genes. Some of the gene do not have a 3'UTR sequence and reported with sequence unavailable but there is FASTA header.
I need to delete or remove all the "sequence unavailable" lines along with their header, i.e. all blank fasta entries.
Please suggest a command with awk or sed command.

Please don't post images of text. Instead, copy/paste the text into your question and use the formatting tools to format it as code. Also, what operating system are you using? We need to know so we can know what tools and what versions of the tools will be available.

You can use awk here, by setting the input (RS) and output (ORS) record separator, this is what defines a "line" in awk, to >. Then, tell it to print records that do not match Sequence unavailable. For instance, using this input file:
$ cat file
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable
>bar1
AGAGATAGAGAGATATA
>asld
Sequence unavailable
>asd
ACTTC

You can do:
$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC
>

That will add a trailing >, so you need one more step to remove it:
$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file | sed 's/>$//'
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC

Alternatively, using the FastaToTbl and TblToFasta scripts I have posted here, you can do:
$ FastaToTbl file | grep -v 'Sequence unavailable' | TblToFasta
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC

Both approaches have the advantage of not needing to read the entire file to memory so can be used on files of any size.

FWIW with GNU awk you could do awk  'BEGIN{RS=">"} {ORS=RT} !/Sequence unavailable/' file instead of piping to sed.

Try this, if you have GNU grep (the default on Linux systems):
grep --no-group-separator -i -B 1 '^[ACGTN]\+$' input.fasta > filtered.fasta

Explanation: find all lines that are composed of only A/C/G/T/N (case insensitive search), and show that line plus the preceding line. Don't put any breaks between matches that are not consecutive.
This should work on both multi-line and single line fasta files.

If there are no line breaks within the sequence data (the format I use) then the following will work.
perl -0777 -nle 'while (/^(>.+\n^[AGCT-]+)\n/gm){print $1,"\n"}' myfile.fa

Test
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable

Output
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA

However If there are line breaks within the sequence data, please use the following ...
perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa

Then pipe to a file,
perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa > myfile_minus_blanks.fa

This works due to a multi-line m
Test ...
>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA
>baz
Sequence unavailable

Out ...
>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA

Note
If you've mixed bases, degenerate nucleotides,
AGCT- should be replaced with AGCTRYWSMKNBDHV (omitting - [not needed])

Note that this will require loading the entire file into memory which might be an issue with large files.

I just processed a 22Gig file with a slightly more complicated command (removing loci). It worked fine.

Thank you very much for your kind response. This is work fro me.

Hi @DrSud great! I use it because its versatile (use it for different stuff). Please do remember to mark an answer "accepted" if you find it helpful and you can also now upvote all the answers given here.

One way using sed:
sed -e '/^>/ { N; }; /Sequence unavailable/ { d; }' file.fa

For each header line, append the next line of input into the pattern space. Then, if the pattern space contains "Sequence unavailable" delete the pattern space and start a new cycle.

Or using awk:
awk '/^>/ { f=$0; next } /^Sequence unavailable/ { next } f { print f; f="" }1' file.fa

For each header line, use a variable to store the current record. We can use next here to force awk to immediately stop processing the current record and move on to the next record, which is just the next line, since, by default, awk's record separator is a newline character. Then, stop processing the current record and move on to the next record if that line begins with "Sequence unavailable". Finally, print the header line if it exists and unset the variable that we used to store it if it does. Once we've printed the header, we can then just print the current record. We can do this by appending 1 to our command string. The 1 on the end just returns true, and by doing so, it performs the default action, which is to print the current record.

There's a bug in the sed code, sed: 1: "/^>/ { N }; /Sequence u ...": extra characters at the end of N command. OSX Darwin

@M__ ah... my mistake. Should be good now I think. Thanks for testing it on OSX

Thanks but same error but as I'm not feeling motivated to test on Linux ...  lets say it an OSX thing, most work on Linux. The awk works fine.

@M__ I don't think I read your error message properly before. Sorry. Hopefully the above will work now..


text:
What does the "filters" argument in "getBM" do?
In the biomaRt package, there is a function getBM which among a few things is useful for mapping between different gene representations. Currently, I'm using it as follows, to map Ensembl transcript IDs to gene names and Ensembl gene IDs.
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       values = as.list(transcripts),
                       mart = mart)

How is the following code different:
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       filters = "ensembl_transcript_id",
                       values = as.list(transcripts),
                       mart = mart)

It's unclear to me from the description of the filters argument here. But just based on results, I get less results when I run the second one.

Hi An Ignorant Wanderer,
Filters allow you to narrow down the dataset to query a specific subset of features within your dataset. A description and tutorial of the filtering steps within BioMart can be found here:
https://www.ensembl.org/info/data/biomart/how_to_use_biomart.html
If you remove filters, you will query the full dataset (all human genes, in the example above). However, this is not recommended for BioMart as large, genome-wide queries will overload the BioMart servers. For queries of this size, we suggest using the Ensembl REST API or FTP site.


text:
Why do I have missing values returned from getBM when converting Ensembl transcript IDs to gene names?
When I run the following R code to convert Ensembl transcript IDs to gene names and Ensembl gene IDs:
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       filters = "ensembl_transcript_id",
                       values = as.list(transcripts),
                       mart = mart)

some of rows returned have NA values for both the gene name and Ensembl gene ID. This is a bit odd because for example, one of the transcript IDs with missing gene IDs is ENST00000414213 but when I google it, I find that it corresponds to gene ENSG00000235143 (according to this). Why am I getting missing results in getBM then? Is it because maybe the transcript is from an old version of the database? If so, is that the only reason or are there other reasons this could be happening?

Hi An Ignorant Wanderer,
Yes, you're right - ENST00000414213 is a deprecated ID from Ensembl release 75.
https://www.ensembl.org/Homo_sapiens/Transcript/Idhistory?t=ENST00000414213
To convert a list of IDs from GRCh37, you will need to use the GRCh37 dataset. You can do this using the following documentation/example:
https://www.ensembl.org/info/data/biomart/biomart_r_package.html


text:
Non-parametric, background-based tests on proteomics data from Proteome Discoverer v2.5?
I am new to proteomics analysis, so apologies if my question is silly! For context, I am working with proteome samples from the postmortem human brain for a case-control study.
Our lab generally receives an output file of Proteome Discoverer* (v2.5) from the proteomics center that we work with. From what I understand, PD perform t-tests on the normalized data, but it is a background-based t-test, so there is a weighting factor is applied that takes into account the distribution of ratios for all other proteins and/or the protein abundances themselves? As such, I am able to see an "Abundance Ratio Adj. P-Value".
My problem is that my data does not meet the assumptions for a parametric test and it has a small sample size, so the t-test is not valid (as far as I can tell, unless there is some sort of log transformation that PD does to the data before the t-test). Therefore, I would like to run a Mann-Whitney U Test. Can this be done in the PD software?
If not, how can I do this background adjustment manually such that I can run a Mann-Whitney U Test but still take into protein abundance like in PD? I do most of my stats in R in case that is relevant.
Thank you so much in advance for your help!
Kelly

Package distributed by ThermoFisher

Overall, this is a good call and the right thing to do. It is worth contact ThermoFisher tech. support to confirm the absence of non-parametric in the packages. Thats what the labs paying for.

Thanks for your reply! I actually did contact Thermo, it doesn't seem as though they have non-parametric options, but I am trying to get information on the specific calculations used for controlling for protein abundance/abundance ratio adjusted p-value. Unfortunately, I am not getting anything helpful thus far :(

Abundance Ratio Adj. P-Value - How is it calculated - the formula and procedure?what do you mean by normalized data?

Small samples size and parametric tests are not mutually exclusive. Specialized software has been developed to leverage sharing information across  measurements to boost power, that is called Empirical Bayes. You can use the limma package for this sort of analysis, contrasting groups. Input is expected to be normalized and on log2-scale. It's an R package hosted at Bioconductor. Don't try a MWU, with small sample size that is massively underpowered.

I appreciate the contribution and yes as a basic principle non-parametric tests have lower statistical power than their parametric alternatives, of course. However a significant result with a non-parametric tests is robust because it is assumption free. As we have seen with DSeq2 recently the problem start when a test is frequently deployed outside its original context, sooner or later someone identifies discrepancies and the whole circus falls over.

Thank you so much for your reply! The reason that I wanted to use a non-parametric test is because the normalized data I receive from PD does not meet the assumptions for a t-test. I am able to do a Mann-Whitney U test, but the problem is that it does not take into account the protein abundance as described above. Do you know if the limma package you proposed can do this?
Thanks again so much for your help.


text:
Only the first file Pair is processed from Channel.fromFilePairs()
This is my first nextflow pipeline and i am struggling with the .fromFilePairs concept of Nextflow.
I want to process paired-end read files and align them against a reference genome. I use Hisat2 so the index files come in form of 8 files with the suffix .ht2. Hisat needs the file path and the filename without ".ht2". All worked well with a single-end version of the process. The reference also has to be passed through a channel because i use a hisat docker-container and otherwise the files are not mounted.
The idea is to pass read-pairs to a preprocessing step and then pass them (as tuples, like a Channel.fromPairs()) into the hisat process and then further downstream.
My problem is all file-pairs are processed during the first step but only the first pair is passed further downstream. All downstream steps finish with the one pair and the rest simply stays in the output-folder of the preprocessing step.
I tried to simply skip the first step to see if the "connection" between the first and second process is not working but the Channel which works for the first step does not for the hisat step. So i think there is a problem with the combination of multiple read-pairs as input but only one reference file?!
process hisat2_paired {
  tag "$sample_id"

  container "${params.container_hisat}"
  publishDir "${params.out_dir}/Nextflow_output/hisat/"

  input:
  tuple val(sample_id), path(input_files)
  path( index_files )

  output:
  path('*.sam') , emit: hisat2_sam_out_ch
  path('*.txt')

  script:
  """
  INDEX=`find -L ${index_files}/ -name "*.1.ht2" | sed 's/.1.ht2//'`
  hisat2 --rna-strandness RF \
          -x \\\$INDEX \
          -1 ${input_files[0]} \
      -2 ${input_files[1]} \
          -S ${sample_id}.sam &> ${sample_id}.hisat_summary.txt
  """

workflow {

Channel
    .fromFilePairs( "${params.reads_folder}/*_{R1,R2}*f*q.gz")
    .ifEmpty { exit 1, "no fastq files found at given path" }
    .view()
    .set{ read_pair_ch }

Channel
    .fromPath( params.hisat2_index )
    .ifEmpty { exit 1, "no hisat2 index files found - path was empty" }
    //.view()
    .set { hisat2_index_ch }

hisat2_paired( read_pair_ch, hisat2_index_ch )

}

I started the workflow with 3 test file-pairs and the Output on the CLI states that only one process is started although if i look at the Channel.view() all 3 file pairs are printed.
executor >  local (1)
[2f/30a4d6] process > hisat2_paired (MC3T3-Ko1_S2_L001) [  0%] 0 of 1
[MC3T3-Ko1_S2_L001, [/mnt/c/Path/to/MC3T3-Ko1_S2_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko1_S2_L001_R2_001.fastq.gz]]
[MC3T3-Ko2_S4_L001, [/mnt/c/Path/to/MC3T3-Ko2_S4_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko2_S4_L001_R2_001.fastq.gz]]
[MC3T3-Ko3_S6_L001, [/mnt/c/Path/to/MC3T3-Ko3_S6_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko3_S6_L001_R2_001.fastq.gz]]

I was already thinking about memory and resources. Since I prototype on my local PC i set the queue size in the config file to 2. I know that my hisat process needs about 7Gb RAM (for Mus.musculus) per process so there should be enough resources. Also the problem was that the local executor is not even trying to start a second process.

So i think there is a problem with the combination of multiple
read-pairs as input but only one reference file?!

Correct. The problem is that the hisat2_index_ch is a queue channel, but what you want is a value channel so that the (file) values can be read an unlimited number of times without being consumed. Assuming params.hisat2_index is intended to supply a glob pattern, you could simply call the collect operator to produce a list of files as the sole emission:
Channel
    .fromPath( params.hisat2_index )
    .collect()
    .set { hisat2_index_ch }

You would then need to update your process definition accordingly. Alternatively, have your workflow create the HISAT2 index files for you. This simplifies things considerably, by requiring only a single reference FASTA file as input. Note that value channels are implicitly created as output for a process whose inputs are all value channels. This means in the example below, hisat2_index.out is also a value channel:
params.publish_dir = './results'

params.container_hisat = 'quay.io/biocontainers/hisat2:2.2.1--h87f3376_4'

params.reads = '*_R{1,2}.fastq.gz'
params.ref_fasta = 'chr1.fasta'

process hisat2_index {

    container params.container_hisat
    cpus 8

    input:
    path fasta

    output:
    tuple val(fasta.baseName), path("${fasta.baseName}.*.ht2")

    """
    hisat2-build \\
        -p ${task.cpus} \\
        "${fasta}" \\
        "${fasta.baseName}"
    """
}

process hisat2_paired {

    tag { sample_id }

    container params.container_hisat
    publishDir "${params.publish_dir}/hisat"

    input:
    tuple val(sample_id), path(input_files)
    tuple val(idx), path('hisat2-idx/*')

    output:
    path('*.sam'), emit: hisat2_sam
    path('*.txt')

    script:
    def (fq1, fq2) = input_files

    """
    hisat2 \\
      --rna-strandness RF \\
      -x "hisat2-idx/${idx}" \\
      -1 ${fq1} \\
      -2 ${fq2} \\
      -S ${sample_id}.sam \\
      &> ${sample_id}.hisat_summary.txt
    """
}

workflow {

    read_pairs = Channel.fromFilePairs( params.reads )
    ref_fasta = file( params.ref_fasta )

    hisat2_index( ref_fasta )

    hisat2_paired( read_pairs, hisat2_index.out )
}

Thank you for the help. I knew it was something simple. Simply adding the .collect() to the hisat_index_ch did the trick. But the process definition needed to be set as path or else the files are not available in the docker container.
By "update the process definition" you meant changing the input from path() to val() right?

@stfn.snk no, sorry, I should have elaborated. Always use the path() qualifier for file input/outputs. The path() qualifier will work for a single file or a collection of files. I just meant that you will need to update how you define your HISAT2 INDEX variable. I don't think the find command will work as expected. It might be better to supply the basename of the index (i.e. a string) using val() in a tuple() like in the example above.

One way, could be to use the map operator to transform the collection, similar to how the find command was intended to work, .e.g Channel.fromPath( params.hisat2_index ).collect().map { tuple( it.first().getBaseName(2), it ) }.set { ... }

I tried to avoid the building process, since it needs considerable computing resources (so i've read) so I used a public available HISAT Index. I think i tried the .baseName function already but it gave me the index filename including the number and the process failed (e.g. genome.1.ht2 became genome.1) So i ended up with the find solution.
I will look further into your suggestions. Thanks again, i really appreciate your help

@stfn.snk So you could use .getBaseName(2) to remove the double extension, or .simpleName if you can be sure the prefix doesn't contain a literal dot/period. Anytime!