text:
How to count RNA sequencereads using custom made python scripts?
I am trying to do RNA seq analysis and my goal is to filter gene counts less than 5 using custom made python scripts. The code chunk goes as follows
#!/usr/bin/python
#Load relevant modules
import sys

#Read in the HTseq gene count output files for each of the 8 samples
inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]
##Determine the genes with sufficient coverage
#Set up a variable to hold genes with sufficient coverage in any of the 8 samples
genelist = []

#For each HTseq file
for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')
    #For each line
    for line in inputfile:
        #Determine the number of reads/transcripts for this gene
        splitline = line.strip().split('\t')
        #Does it pass the threshold?
        if int(splitline[1]) >= 5:
            #If so, add it to list
            if not splitline[0] in genelist:
                genelist.append(splitline[0])

    #Close the file
    inputfile.close()

#Write out the list of sufficiently expressed genes
outputlist = open(sys.argv[9], 'w')

for gene in genelist:
    outputlist.write(gene + '\n')

outputlist.close()
##Filter each of the HTseq files for only sufficiently expressed genes across our samples
#For each HTseq file

for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')

#Open an output 'filtered' file
outputfile = open(sys.argv[i+10], 'w')

#For each line
for line in inputfile:
    #Determine the gene
    splitline = line.strip().split('\t')

    #Is it in our list?
    if splitline[0] in genelist:
        #If so, write out the original line
        outputfile.write(line)

#Close the files
inputfile.close()
outputfile.close()

python filtercounts.py 'wth1count.txt' 'wth2count.txt' 'wth3count.txt' 'wth4count.txt' '99h1count.txt' '99h2count.txt' '99h3count.txt' '99h4count.txt' 'genelist.txt' 'wth1filtered.txt' 'wth2filtered.txt' 'wth3filtered.txt' 'wth4filtered.txt'  '99h1filtered.txt' '99h2filtered.txt' '99h3filtered.txt' '99h4filtered.txt'
The program works fine but at the end of the code execution I am not getting the last two output files '99h3filtered.txt' '99h4filtered.txt'.

If you want to provide a small sample data I will show you hwo to do it in pandas

What software are you using for the next step?  DESeq in R will do this filtering for you in one line.
I think a more typical way to do this would be to make your entire gene count matrix, and then remove rows where the max value for that row is < 5.  Your way, you have to go through each file twice; once to know which genes to keep, and then again to actually remove lines.


text:
How to annotate optimally a fungal genome without RNA-seq evidence?
Genome information：

~50M nt
2300+ contigs
No pre-trained parameters in Augustus
There are several well annotated RefSeq genomes of this genus.
There are several RNA-seq data of this genus, but not in this species.

My current strategy:
EST evidence:

I selected many RNA-seq data that include as many species as
possible, and de novo assembled these data using Trinity.
Only the longest transcripts were saved.
All transcript files were catenated into a Total_transcript.fasta file, and cd-hit-est was used to reduce data size with options -c 0.8 -n 5 -s 0.8.

Protein evidence

All fungal proteins from SwissProt was download and reduced using cd-hit with options -c 0.8 -n 5 -s 0.8.

Proteinsets from RefSeq genomes of this genus were download, catenated, reduced using same methods.

I planed to annotate genome using Maker2 with above Evidence.
I have a doubt about my strategy: Maker2 calls Augustus to ab-initio annotation, but no pre-trained model. How should I do to solve this?
Anyone else could share your valuable suggestions? Thanks anyway.

The Maker documentation does include some instructions for training ab initio gene predictors, but it assumes an abundant EST database is available. (Assumptions about what kinds of sequences will be available for a draft genome assembly have changed drastically in the last decade.)
It may be worth exploring an iterative approach in any case. You can do a first-pass annotation with Maker using the evidence you discussed, along with ab initio predictions from Augustus (and I'd recommend SNAP as well). Each gene model will be scored by Maker with the annotation edit distance (AED), representing how well the gene model agrees with the evidence (0.0 is best). From this first pass, if you can identify a few hundred reliable gene models with good AED scores (or maybe even as many as one or two thousand), that should be plenty of training data. You could use these gene models to train a new model for Augustus, using Maker's instructions or Augustus' own instruction. (In my experience, it was quite a bit of work in either case.)

source: Stack Exchange

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

sort: votes


text:
FastqCleaner in R, not recognizing paired-end fastq file?
I've recently downloaded the FastqCLeaner to process paired-end fastq files locally using R (4.0.3). I'm having trouble with the shiny app. It is not uploading my second fastq file, and is repeated asking for file 1. Can someone help me? is there another package i can use for trimming and alignment in R? I've tried GREP2 but had massive installation issues, i wasted 2 days trying to install in R.

Ditch R, it's not the proper environment for handling reads. R is for statistical analysis and plots and absolutely nothing else.

Can you clarify what tasks need to be done on those files? There is most likely a simple command line tool that handles this in like no time without the need to bother with R-ish tools here.

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


text:
How to write multi-sequence alignment data frame as a fasta format fast
I have a data frame where the row names are the species names, and each column can either be an Amino-Acid character or "-".
I wish to write it in a fasta format. Simple example of the format is:
>UniRef90_A0A1V4QIJ3
------------------------------------------------------------
--------------------------MDISTFSP-----PVAG-----------VE----
------ITAF-------------V------------------------------------
LN--------G-------------------------------------------------
----EKIPA-----------------------------------------IVLSKD----
---LNPFL------------------Q-ELLDQEPAC-----------------------
------------EHDIGTCYG------------------------DA-------------
--------------------------------------------------PC-----IS-
-----------------------GPDAWNHSGFQC----Q---------VC---------
------------GLVKSWNL-------------------VGDFVVYRQE-----------
>UniRef90_A0A0L9V5F8
------------------------------------------------------------
RFHTLFRNEYGHLRVLQRFDQRSKQIQNLENYRLVEFKSKPNTLLLPHHADADFLLVVLN
GRALLTLVNP-------------D---------------------------------G--
RDSYILEQGHA-------------------------------------------------
----QKIPAGTIFFLVNPNDNENLRI------------------------IKIATPINNP
HRFQDFFL-------------SSTEAQQSYLQG--FS-----------------------
----------------KNVLEASFDSEFKEINRVLFGEEGQQQQGEESQQEGVIVELERE
QIRELIKHA-----------------------------------------KSSSRRSLS-
-----------------------SQ-----DEPFNLRNRKP--------IY--SN-----
------------KF-GRWYEITPEKNP------------Q--------------------

I can do it using the following for loop:
for (n in row.names(df)){
    ss = do.call(paste, c(as.list(mu_df[1, ]), sep = ""))
    write.fasta(ss,
                  names=n,
                  as.string=TRUE,
                  open="a",
                  new_file_path)
}

Is there a faster way of writing it? I have to write millions of those files and I wish to optimize it.
In order to create the data frame:
structure(list(V1 = c("-", "-"), V2 = c("-", "-"), V3 = c("-",
"-"), V4 = c("-", "-"), V5 = c("-", "-"), V6 = c("-", "-"), V7 = c("-",
"-"), V8 = c("-", "-"), V9 = c("-", "-"), V10 = c("-", "-"),
    V11 = c("-", "-"), V12 = c("-", "-"), V13 = c("-", "-"),
    V14 = c("-", "-"), V15 = c("-", "-"), V16 = c("-", "-"),
    V17 = c("-", "-"), V18 = c("-", "-"), V19 = c("-", "-"),
    V20 = c("-", "-"), V21 = c("-", "-"), V22 = c("-", "-"),
    V23 = c("-", "-"), V24 = c("-", "-"), V25 = c("-", "-"),
    V26 = c("-", "-"), V27 = c("-", "-"), V28 = c("-", "-"),
    V29 = c("-", "-"), V30 = c("-", "-"), V31 = c("-", "-"),
    V32 = c("-", "-"), V33 = c("-", "-"), V34 = c("-", "-"),
    V35 = c("-", "-"), V36 = c("-", "-"), V37 = c("-", "-"),
    V38 = c("-", "-"), V39 = c("-", "-"), V40 = c("-", "-"),
    V41 = c("-", "-"), V42 = c("-", "-"), V43 = c("-", "-"),
    V44 = c("-", "-"), V45 = c("-", "-"), V46 = c("-", "-"),
    V47 = c("-", "-"), V48 = c("-", "-"), V49 = c("-", "-"),
    V50 = c("-", "-"), V51 = c("-", "-"), V52 = c("-", "-"),
    V53 = c("-", "-"), V54 = c("-", "-"), V55 = c("-", "-"),
    V56 = c("-", "-"), V57 = c("-", "-"), V58 = c("-", "-"),
    V59 = c("-", "-"), V60 = c("-", "-"), V61 = c("-", "R"),
    V62 = c("-", "F"), V63 = c("-", "H"), V64 = c("-", "T"),
    V65 = c("-", "L"), V66 = c("-", "F"), V67 = c("-", "R"),
    V68 = c("-", "N"), V69 = c("-", "E"), V70 = c("-", "Y"),
    V71 = c("-", "G"), V72 = c("-", "H"), V73 = c("-", "L"),
    V74 = c("-", "R"), V75 = c("-", "V"), V76 = c("-", "L"),
    V77 = c("-", "Q"), V78 = c("-", "R"), V79 = c("-", "F"),
    V80 = c("-", "D"), V81 = c("-", "Q"), V82 = c("-", "R"),
    V83 = c("-", "S"), V84 = c("-", "K"), V85 = c("-", "Q"),
    V86 = c("-", "I"), V87 = c("M", "Q"), V88 = c("D", "N"),
    V89 = c("I", "L"), V90 = c("S", "E"), V91 = c("T", "N"),
    V92 = c("F", "Y"), V93 = c("S", "R"), V94 = c("P", "L"),
    V95 = c("-", "V"), V96 = c("-", "E"), V97 = c("-", "F"),
    V98 = c("-", "K"), V99 = c("-", "S"), V100 = c("P", "K"),
    V101 = c("V", "P"), V102 = c("A", "N"), V103 = c("G", "T"
    ), V104 = c("-", "L"), V105 = c("-", "L"), V106 = c("-",
    "L"), V107 = c("-", "P"), V108 = c("-", "H"), V109 = c("-",
    "H"), V110 = c("-", "A"), V111 = c("-", "D"), V112 = c("-",
    "A"), V113 = c("-", "D"), V114 = c("-", "F"), V115 = c("V",
    "L"), V116 = c("E", "L"), V117 = c("-", "V"), V118 = c("-",
    "V"), V119 = c("-", "L"), V120 = c("-", "N"), V121 = c("-",
    "G"), V122 = c("-", "R"), V123 = c("-", "A"), V124 = c("-",
    "L"), V125 = c("-", "L"), V126 = c("-", "T"), V127 = c("I",
    "L"), V128 = c("T", "V"), V129 = c("A", "N"), V130 = c("F",
    "P"), V131 = c("-", "-"), V132 = c("-", "-"), V133 = c("-",
    "-"), V134 = c("-", "-"), V135 = c("-", "-"), V136 = c("-",
    "-"), V137 = c("-", "-"), V138 = c("-", "-"), V139 = c("-",
    "-"), V140 = c("-", "-"), V141 = c("-", "-"), V142 = c("-",
    "-"), V143 = c("-", "-"), V144 = c("V", "D"), V145 = c("-",
    "-"), V146 = c("-", "-"), V147 = c("-", "-"), V148 = c("-",
    "-"), V149 = c("-", "-"), V150 = c("-", "-"), V151 = c("-",
    "-"), V152 = c("-", "-"), V153 = c("-", "-"), V154 = c("-",
    "-"), V155 = c("-", "-"), V156 = c("-", "-"), V157 = c("-",
    "-"), V158 = c("-", "-"), V159 = c("-", "-"), V160 = c("-",
    "-"), V161 = c("-", "-"), V162 = c("-", "-"), V163 = c("-",
    "-"), V164 = c("-", "-"), V165 = c("-", "-"), V166 = c("-",
    "-"), V167 = c("-", "-"), V168 = c("-", "-"), V169 = c("-",
    "-"), V170 = c("-", "-"), V171 = c("-", "-"), V172 = c("-",
    "-"), V173 = c("-", "-"), V174 = c("-", "-"), V175 = c("-",
    "-"), V176 = c("-", "-"), V177 = c("-", "-"), V178 = c("-",
    "G"), V179 = c("-", "-"), V180 = c("-", "-"), V181 = c("L",
    "R"), V182 = c("N", "D"), V183 = c("-", "S"), V184 = c("-",
    "Y"), V185 = c("-", "I"), V186 = c("-", "L"), V187 = c("-",
    "E"), V188 = c("-", "Q"), V189 = c("-", "G"), V190 = c("-",
    "H"), V191 = c("G", "A"), V192 = c("-", "-"), V193 = c("-",
    "-"), V194 = c("-", "-"), V195 = c("-", "-"), V196 = c("-",
    "-"), V197 = c("-", "-"), V198 = c("-", "-"), V199 = c("-",
    "-"), V200 = c("-", "-"), V201 = c("-", "-"), V202 = c("-",
    "-"), V203 = c("-", "-"), V204 = c("-", "-"), V205 = c("-",
    "-"), V206 = c("-", "-"), V207 = c("-", "-"), V208 = c("-",
    "-"), V209 = c("-", "-"), V210 = c("-", "-"), V211 = c("-",
    "-"), V212 = c("-", "-"), V213 = c("-", "-"), V214 = c("-",
    "-"), V215 = c("-", "-"), V216 = c("-", "-"), V217 = c("-",
    "-"), V218 = c("-", "-"), V219 = c("-", "-"), V220 = c("-",
    "-"), V221 = c("-", "-"), V222 = c("-", "-"), V223 = c("-",
    "-"), V224 = c("-", "-"), V225 = c("-", "-"), V226 = c("-",
    "-"), V227 = c("-", "-"), V228 = c("-", "-"), V229 = c("-",
    "-"), V230 = c("-", "-"), V231 = c("-", "-"), V232 = c("-",
    "-"), V233 = c("-", "-"), V234 = c("-", "-"), V235 = c("-",
    "-"), V236 = c("-", "-"), V237 = c("-", "-"), V238 = c("-",
    "-"), V239 = c("-", "-"), V240 = c("-", "-"), V241 = c("-",
    "-"), V242 = c("-", "-"), V243 = c("-", "-"), V244 = c("-",
    "-"), V245 = c("E", "Q"), V246 = c("K", "K"), V247 = c("I",
    "I"), V248 = c("P", "P"), V249 = c("A", "A"), V250 = c("-",
    "G"), V251 = c("-", "T"), V252 = c("-", "I"), V253 = c("-",
    "F"), V254 = c("-", "F"), V255 = c("-", "L"), V256 = c("-",
    "V"), V257 = c("-", "N"), V258 = c("-", "P"), V259 = c("-",
    "N"), V260 = c("-", "D"), V261 = c("-", "N"), V262 = c("-",
    "E"), V263 = c("-", "N"), V264 = c("-", "L"), V265 = c("-",
    "R"), V266 = c("-", "I"), V267 = c("-", "-"), V268 = c("-",
    "-"), V269 = c("-", "-"), V270 = c("-", "-"), V271 = c("-",
    "-"), V272 = c("-", "-"), V273 = c("-", "-"), V274 = c("-",
    "-"), V275 = c("-", "-"), V276 = c("-", "-"), V277 = c("-",
    "-"), V278 = c("-", "-"), V279 = c("-", "-"), V280 = c("-",
    "-"), V281 = c("-", "-"), V282 = c("-", "-"), V283 = c("-",
    "-"), V284 = c("-", "-"), V285 = c("-", "-"), V286 = c("-",
    "-"), V287 = c("-", "-"), V288 = c("-", "-"), V289 = c("-",
    "-"), V290 = c("-", "-"), V291 = c("I", "I"), V292 = c("V",
    "K"), V293 = c("L", "I"), V294 = c("S", "A"), V295 = c("K",
    "T"), V296 = c("D", "P"), V297 = c("-", "I"), V298 = c("-",
    "N"), V299 = c("-", "N"), V300 = c("-", "P"), V301 = c("-",
    "H"), V302 = c("-", "R"), V303 = c("-", "F"), V304 = c("L",
    "Q"), V305 = c("N", "D"), V306 = c("P", "F"), V307 = c("F",
    "F"), V308 = c("L", "L"), V309 = c("-", "-"), V310 = c("-",
    "-"), V311 = c("-", "-"), V312 = c("-", "-"), V313 = c("-",
    "-"), V314 = c("-", "-"), V315 = c("-", "-"), V316 = c("-",
    "-"), V317 = c("-", "-"), V318 = c("-", "-"), V319 = c("-",
    "-"), V320 = c("-", "-"), V321 = c("-", "-"), V322 = c("-",
    "S"), V323 = c("-", "S"), V324 = c("-", "T"), V325 = c("-",
    "E"), V326 = c("-", "A"), V327 = c("Q", "Q"), V328 = c("-",
    "Q"), V329 = c("E", "S"), V330 = c("L", "Y"), V331 = c("L",
    "L"), V332 = c("D", "Q"), V333 = c("Q", "G"), V334 = c("E",
    "-"), V335 = c("P", "-"), V336 = c("A", "F"), V337 = c("C",
    "S"), V338 = c("-", "-"), V339 = c("-", "-"), V340 = c("-",
    "-"), V341 = c("-", "-"), V342 = c("-", "-"), V343 = c("-",
    "-"), V344 = c("-", "-"), V345 = c("-", "-"), V346 = c("-",
    "-"), V347 = c("-", "-"), V348 = c("-", "-"), V349 = c("-",
    "-"), V350 = c("-", "-"), V351 = c("-", "-"), V352 = c("-",
    "-"), V353 = c("-", "-"), V354 = c("-", "-"), V355 = c("-",
    "-"), V356 = c("-", "-"), V357 = c("-", "-"), V358 = c("-",
    "-"), V359 = c("-", "-"), V360 = c("-", "-"), V361 = c("-",
    "-"), V362 = c("-", "-"), V363 = c("-", "-"), V364 = c("-",
    "-"), V365 = c("-", "-"), V366 = c("-", "-"), V367 = c("-",
    "-"), V368 = c("-", "-"), V369 = c("-", "-"), V370 = c("-",
    "-"), V371 = c("-", "-"), V372 = c("-", "-"), V373 = c("E",
    "-"), V374 = c("H", "-"), V375 = c("D", "-"), V376 = c("I",
    "-"), V377 = c("G", "K"), V378 = c("T", "N"), V379 = c("C",
    "V"), V380 = c("Y", "L"), V381 = c("G", "E"), V382 = c("-",
    "A"), V383 = c("-", "S"), V384 = c("-", "F"), V385 = c("-",
    "D"), V386 = c("-", "S"), V387 = c("-", "E"), V388 = c("-",
    "F"), V389 = c("-", "K"), V390 = c("-", "E"), V391 = c("-",
    "I"), V392 = c("-", "N"), V393 = c("-", "R"), V394 = c("-",
    "V"), V395 = c("-", "L"), V396 = c("-", "F"), V397 = c("-",
    "G"), V398 = c("-", "E"), V399 = c("-", "E"), V400 = c("-",
    "G"), V401 = c("-", "Q"), V402 = c("-", "Q"), V403 = c("-",
    "Q"), V404 = c("-", "Q"), V405 = c("-", "G"), V406 = c("D",
    "E"), V407 = c("A", "E"), V408 = c("-", "S"), V409 = c("-",
    "Q"), V410 = c("-", "Q"), V411 = c("-", "E"), V412 = c("-",
    "G"), V413 = c("-", "V"), V414 = c("-", "I"), V415 = c("-",
    "V"), V416 = c("-", "E"), V417 = c("-", "L"), V418 = c("-",
    "E"), V419 = c("-", "R"), V420 = c("-", "E"), V421 = c("-",
    "Q"), V422 = c("-", "I"), V423 = c("-", "R"), V424 = c("-",
    "E"), V425 = c("-", "L"), V426 = c("-", "I"), V427 = c("-",
    "K"), V428 = c("-", "H"), V429 = c("-", "A"), V430 = c("-",
    "-"), V431 = c("-", "-"), V432 = c("-", "-"), V433 = c("-",
    "-"), V434 = c("-", "-"), V435 = c("-", "-"), V436 = c("-",
    "-"), V437 = c("-", "-"), V438 = c("-", "-"), V439 = c("-",
    "-"), V440 = c("-", "-"), V441 = c("-", "-"), V442 = c("-",
    "-"), V443 = c("-", "-"), V444 = c("-", "-"), V445 = c("-",
    "-"), V446 = c("-", "-"), V447 = c("-", "-"), V448 = c("-",
    "-"), V449 = c("-", "-"), V450 = c("-", "-"), V451 = c("-",
    "-"), V452 = c("-", "-"), V453 = c("-", "-"), V454 = c("-",
    "-"), V455 = c("-", "-"), V456 = c("-", "-"), V457 = c("-",
    "-"), V458 = c("-", "-"), V459 = c("-", "-"), V460 = c("-",
    "-"), V461 = c("-", "-"), V462 = c("-", "-"), V463 = c("-",
    "-"), V464 = c("-", "-"), V465 = c("-", "-"), V466 = c("-",
    "-"), V467 = c("-", "-"), V468 = c("-", "-"), V469 = c("-",
    "-"), V470 = c("-", "-"), V471 = c("P", "K"), V472 = c("C",
    "S"), V473 = c("-", "S"), V474 = c("-", "S"), V475 = c("-",
    "R"), V476 = c("-", "R"), V477 = c("-", "S"), V478 = c("I",
    "L"), V479 = c("S", "S"), V480 = c("-", "-"), V481 = c("-",
    "-"), V482 = c("-", "-"), V483 = c("-", "-"), V484 = c("-",
    "-"), V485 = c("-", "-"), V486 = c("-", "-"), V487 = c("-",
    "-"), V488 = c("-", "-"), V489 = c("-", "-"), V490 = c("-",
    "-"), V491 = c("-", "-"), V492 = c("-", "-"), V493 = c("-",
    "-"), V494 = c("-", "-"), V495 = c("-", "-"), V496 = c("-",
    "-"), V497 = c("-", "-"), V498 = c("-", "-"), V499 = c("-",
    "-"), V500 = c("-", "-"), V501 = c("-", "-"), V502 = c("-",
    "-"), V503 = c("-", "-"), V504 = c("G", "S"), V505 = c("P",
    "Q"), V506 = c("D", "-"), V507 = c("A", "-"), V508 = c("W",
    "-"), V509 = c("N", "-"), V510 = c("H", "-"), V511 = c("S",
    "D"), V512 = c("G", "E"), V513 = c("F", "P"), V514 = c("Q",
    "F"), V515 = c("C", "N"), V516 = c("-", "L"), V517 = c("-",
    "R"), V518 = c("-", "N"), V519 = c("-", "R"), V520 = c("Q",
    "K"), V521 = c("-", "P"), V522 = c("-", "-"), V523 = c("-",
    "-"), V524 = c("-", "-"), V525 = c("-", "-"), V526 = c("-",
    "-"), V527 = c("-", "-"), V528 = c("-", "-"), V529 = c("-",
    "-"), V530 = c("V", "I"), V531 = c("C", "Y"), V532 = c("-",
    "-"), V533 = c("-", "-"), V534 = c("-", "S"), V535 = c("-",
    "N"), V536 = c("-", "-"), V537 = c("-", "-"), V538 = c("-",
    "-"), V539 = c("-", "-"), V540 = c("-", "-"), V541 = c("-",
    "-"), V542 = c("-", "-"), V543 = c("-", "-"), V544 = c("-",
    "-"), V545 = c("-", "-"), V546 = c("-", "-"), V547 = c("-",
    "-"), V548 = c("-", "-"), V549 = c("-", "-"), V550 = c("-",
    "-"), V551 = c("-", "-"), V552 = c("-", "-"), V553 = c("G",
    "K"), V554 = c("L", "F"), V555 = c("V", "-"), V556 = c("K",
    "G"), V557 = c("S", "R"), V558 = c("W", "W"), V559 = c("N",
    "Y"), V560 = c("L", "E"), V561 = c("-", "I"), V562 = c("-",
    "T"), V563 = c("-", "P"), V564 = c("-", "E"), V565 = c("-",
    "K"), V566 = c("-", "N"), V567 = c("-", "P"), V568 = c("-",
    "-"), V569 = c("-", "-"), V570 = c("-", "-"), V571 = c("-",
    "-"), V572 = c("-", "-"), V573 = c("-", "-"), V574 = c("-",
    "-"), V575 = c("-", "-"), V576 = c("-", "-"), V577 = c("-",
    "-"), V578 = c("-", "-"), V579 = c("-", "-"), V580 = c("V",
    "Q"), V581 = c("G", "-"), V582 = c("D", "-"), V583 = c("F",
    "-"), V584 = c("V", "-"), V585 = c("V", "-"), V586 = c("Y",
    "-"), V587 = c("R", "-"), V588 = c("Q", "-"), V589 = c("E",
    "-"), V590 = c("-", "-"), V591 = c("-", "-"), V592 = c("-",
    "-"), V593 = c("-", "-"), V594 = c("-", "-"), V595 = c("-",
    "-"), V596 = c("-", "-"), V597 = c("-", "-"), V598 = c("-",
    "-"), V599 = c("-", "-"), V600 = c("-", "-")), class = "data.frame", row.names = c("UniRef90_A0A1V4QIJ3",
"UniRef90_A0A0L9V5F8"))

Update
It seems that part of the desired behavior is to write the sequences in the alignment separately. I am not sure that I understand why it is necessary to write them separately. Generally it is inadvisable to write many small files rather than one big file, based purely on issues of performance; you are likely to hit the I/O bound on performance very quickly, and that will be much slower.
One solution that has been suggested previously for this kind of problem is to write the whole alignment to a single file, and then post-process that filed in a second step to extract the sequences of interest. This could be as trivial as a samtools faidx call from the shell, or you could work with it in seqinr as suggested in the previous answer.
Of course, it seems like it would be simpler for most purposes to just have the big alignment file and use e.g. samtools indexing / faidx / biopython / whatever to extract the sequences of interest. But it seems that that is specifically not desired in this case.
Original answer
I would strongly suggest using the seqinr package, which is optimized for solving exactly these problems. See e.g. this similar question.
You will be able to read in an alignment (or whatever), cast it to as.matrix.alignment, do whatever your analysis is, then cast it back to as.alignment, which can then be passed to write.fasta to write the files.
Whenever possible it is preferable to work with such libraries, which handle these utility functions well, fast, and easy.

Agree, but because I wish to write a lot smaller fasta based on the original one, I was wondering if there might be faster way.

@DavidS I see, it had not occurred to me that anyone would want many small "aligned" fasta files. I have updated the answer with one suggestion, but I think that this is just going to be slow based on file I/O no matter what if you want to write large numbers of files. Possibly someone has something clever though.

So I start from a big MSA, but I wish to do some sort of "sub-sampling" of it to create a smaller files which will still holds valuable information to use further in a deep learning model.

@DavidS Have you considered directly subsampling from the alignment as an in-memory data structure? I think that will definitely be faster than writing files and then reading them in again for the model. If the issue is access to the data structure across different processes, you could write the full alignment as an RData file using save(), and then load() it into each process. (assuming the DL model is in R; if it's in e.g. Python, I am sure that analogous alignment handling is possible there, but that's a different question.)

I considered it, but the problem is that by just sampling from the MSA your might sample species that are too similar to the reference. So I apply some logic to select k samples which are farther apart by some criteria.

@DavidS So why not put that logic into an in-memory sampling regime? Granted, we're now getting into an unrelated implementation discussion, which is mostly my fault. I'll just conclude by reiterating 1) I/O of many files intrinsically makes this (relatively) slow, 2) it is not clear why I/O of many files is necessary, from the high-level discussion we're having.


text:
Required input files for StringTie2
I am working on a virtual project for WGS combined with RNA seq for annotation. The RNA will be sequenced using PacBio Isoseq (Sequel II, HiFi-reads). With some research, I found that StringTie2 is useful for this. In their release paper , they sometimes mention previous alignment or mapping with TopHat or minihat, but I don't find a manual where they say the same. Now I don't know what is supposed to be the input of StringTie2: raw RNA reads, or a mapped BAM-file?
Thank you!

What exactly is unclear after reading the manual: https://github.com/gpertea/stringtie It is mapped reads that are input. I do not know what minihat is, you probably mean minimap2, which is probably a good choice for long reads. For sure don't use tophat. It is old and deprecated, for short reads use something like STAR or hisat2.

I was not sure if StringTie2 has the same requirements, but ok, it's clear that I should use mapped reads. I indeed meant minimap2, must've confused it with tophat ;) Thank you!


text:
How to search for high coverage SRA entries
The Question
I want to find high coverage SRA entries, e.g., above 100x.
I guess the best way is to use https://www.ncbi.nlm.nih.gov/sra with an appropriate search term. I don't mind if the search results contain some "false-positives". (I.e., if some of the found entries have coverage above 80x, while I aimed for 100x, that's fine.)
I want to avoid using Entrez.esummary on a very large number of SRA entries (e.g., through Biopython), as it would be significantly slower (if I understand correctly).
How can I do that?

My Attempt
I tried using the Mbases field.Assuming the reference genome is at most 100kbp (e.g., in case I am interested in viruses), I thought a good approximation would be to search for SRA entries with more than 100*100k bases, i.e., at least 10M bases.
According to https://www.ncbi.nlm.nih.gov/books/NBK3837/, I thought I could use the search term 10:9999999[Mbases], but disappointingly, this didn't work. Just to make sure my syntax is right, I verified that 100:101[ReadLength] works. So it seems that Mbases doesn't support the range syntax.

If processing a huge table with lots of "false positives" is OK, then you could enter your search term, from top right choose "Send to" > "File" Format: "RunInfo".


text:
Is there a baby Hello World example using BIGstack?
Intel and the Broad Institute have developed BIGstack, described in this Broad-Intel Genomics stack whitepaper.
Eventually I would love to perform a variant calling on BIGstack, but it isn't clear how to run a “simple” example of BIGstack on Xeon with Altera FPGA.
Since I have my own old Altera DE2 FPGA board and I would like to run the simplest possible toy example which is possible to run on my old Altera DE2 board. Is it possible to use my Altera DE2 to run something simple on top of BIGstack?

Not sure if this helps but if you are going for FPGA analysis, have you checked Nvidia Parabricks Pipelines? Maybe your FPGA board would work with it.


text:
Can't locate object method "new" via package Bio::Tools::Primer3

Hi there I am still a novice at this. I'm trying to run this perl script called FullSSR https://sourceforge.net/projects/fullssr, by following the instructions in their 'README.TXT' https://sourceforge.net/projects/fullssr/files/.
I did everything as instructed by the app developer(s) which included among others, replacing the Primer3.pm file provided in the *.tar.gz file into /usr/share/perl5/Bio/Tools/Run/; opening the config.txt and specifying the path of primer3_core (file) and primer3_config/ (folder); placing the sequence file I wanted to analyse into the FullSSR folder and executing the perl script using the command line: perl 'FullSSR\(ver-1.1\).p; -i sequence.fasta -m 5.
However, it was to no avail. Instead, I got the message "Can't locate object method "new" via package Bio::Tools::Primer3" (perhaps you forgot to load "Bio::Tools::Primer3?) at /usr/share/perl5/Bio/Tools/Run/Primer3.pm line 497.
I am going to try to install/load Bio::Tools::Primer3 with the cpan module (cpanm) and will see if that does it.
Perhaps someone out there could suggest a helpful (easy to understand) manual on CPAN?!

I managed to figure out the answer to my question. I installed Bio::Tools::Primer3 by typing the following into terminal: 'sudo cpanm Bio::Tools::Primer3 --force' to install it. It only works when I used 'sudo' and 'force'

I managed to figure out the answer to my question. I installed Bio::Tools::Primer3 by typing the following into terminal: 'sudo cpanm Bio::Tools::Primer3 --force' to install it. It only works when I used 'sudo' and 'force'


text:
Nvidia Parabrick fq2bam pipeline error - No such file or directory
I am trying to implement Nvidia-parabrick fq2bam pipeline to process my WES dataset. The reason for opting this pipeline is the huge number of samples (1004) which take many months to process for somatic mutation identification. THe GPU accelerated pipeline will be more helpful to get somatic mutations in a quick manner. After pulling parabrick docker image from here, the following command works fine and return the options :

docker run -it nvcr.io/hpc/parabricks:v2.5.0 -h

But whenever I give the following command to run fq2bam pipeline, I am getting an error of No such file or directory:

docker run -it nvcr.io/hpc/parabricks:v2.5.0 fq2bam --ref ucsc.hg19.fasta --in-fq sample_1.fastq.gz sample_2.fastq.gz --knownSites 1000G_phase1.indels.hg19.sites.vcf.gz --out-bam sample_dups_gpu.bam out-recal-file sample_recal_gpu.txt --tmp-dir /tmp

The error is as follows:
Traceback (most recent call last):
File "/parabricks/run_pipeline.py", line 7, in <module>
sys.exit(PB.pb_main())
File "PB.pyx", line 846, in PB.pb_main
File "/parabricks/pbargs.py", line 1287, in getArgs
return PBRun(sys.argv)
File "/parabricks/pbargs.py", line 392, in __init__
self.runArgs = getattr(self, args.command)(argList)
File "/parabricks/pbargs.py", line 447, in fq2bam
self.addToParser(bsmb_parser_sysgroup, sysOptionGenerator().allOptions)
File "/parabricks/pbargs.py", line 50, in __init__
PBOption(category="sysOption", name="--num-gpus", default=GetNumGPUs(), typeName=int, helpStr="Number of GPUs to use for a run")
File "/parabricks/pbutils.py", line 48, in GetNumGPUs
output = subprocess.check_output(["nvidia-smi", "-L"], universal_newlines=True)
File "/usr/lib/python2.7/subprocess.py", line 567, in check_output
process = Popen(stdout=PIPE, *popenargs, **kwargs)
File "/usr/lib/python2.7/subprocess.py", line 711, in __init__
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1343, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory

As I am using the docker image as it is,it is not clear to me which file is missing. Kindly suggest. Thanks.

Just a tip - if you edit your title to make it more informative (e.g. Nvidia Parabrick error - No such file or directory) it helps people answer the questoin and also helps other people who have the same error.

Thanks. I have approved the edit in question title to make it more informative.

My query was answered by NVIDIA Parabrick developer team. The installation instruction are available at https://forums.developer.nvidia.com/uploads/short-url/t5uXsZgdMXx0G2wDcYamFl9sO2x.pdf.


text:
Clarification regarding (deepTools) bamCoverage output
Im using bamCoverage on a file and ive set the bin size to 1 but also normalizing using RPKM for defined regions using a bedgraph file.  The output from bamCoverage however is in binned segments each with a (RPKM) score.  Am I correct in interpreting then that
1 - the output bedgraph tells me the (total) RPKM score for that particular bin?
2 - to get the total RPKM score for the region defined in the bedgraph I need to sum all the RPKM scores given for that particular region?

If I'm not mistaken, RPKM includes a "normalization" by unit of length, so it does not make sense to sum these values along a region. A weighted average makes more sense (were the weights would be the bin widths).

Yes
The total RPKM for the region in the bedgraph file is the value already present in the bedgraph file. If you want to estimate the RPKM for a larger region composed of multiple bins then you can average the per-bin RPKM values. Note that adjacent bins with the same value are merged by deepTools, so you will need to account for this.

Note that RPKM is NOT a robust normalization and that I should remove it from deepTools so people stop using it.

@OP If you want a more robust normalization I suggest (disclaimer, I am the author of the post) https://www.biostars.org/p/413626/#414440


text:
Can anyone point me towards a HiC Dataset with domains of high confidence?
is there a dataset online containing a HiC contact matrix and file with TAD boundary locations, which were generated, either by hand or automated, with high accuracy? Although I already used programs such as SpectralTAD or GMAP to generate TADs myself, a dataset which was computed with higher accuracy would be nice to check my project against.
Preferred attributes would be:

human genome
whole genome
10 kb resolution or higher

Thank you very much in advance!
edit: One dataset I found is http://chromosome.sdsc.edu/mouse/hi-c/download.html, which includes multiple 10 kb datasets with TADs called.

I guess your best chance is one of the in-situ HiC papers from the Lieberman-Aiden lab. There are not too many HiC datasets with 10kb-or-better resolution out there for human. I would check this one:
A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping


text:
Open chromatin region : nucleosome , mono-nucleosome , di-nucleosome
I m following this tutorial for my ATAC seq analysis. Here in of the section they show profile of open chromatin region plots.
**Plotting insert sizes with Greenleaf open, mono- and di-nucleosome profiles
This looks very similar to the image from the Greenleaf paper.

We can now annotate our nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) as in the Greenleaf study.**

Now Since i have done the differential accessibility analysis. So from that information how do i find out region which are "nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) ".
My objective

I would like to see this "nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) " from my differential accessible region which was done from readcount and then I ran deseq2.
My goal is not to plot the way they have done rather find those regions mentioned above and see it between different conditions.

Now based on the tutorial of that section where they do for each bam file.
Do i need to run the analysis on each bam file or each condition or Can I find out the above from my deseq2 analysis which has regions and those regions I have annotated it .

It's impossible to find nucleosome free and such regions without a BAM file. Your only options are:

Filter the BAM files to contain only the size-range of interest and rerun the counting and statistics.
Find or write a counting program that can track fragment sizes.

BTW, those size ranges could use some adjustments, you're throwing away all fragments between 100 and 180 bases using their ranges.

Thank you devon . "BTW, those size ranges could use some adjustments, you're throwing away all fragments between 100 and 180 bases using their ranges." how to do it add some region?

Looks at a size distribution plot and adjust accordingly.


text:
What's the meaning of Txn factor?
The transcription factor track in UCSC is named as Txn factor. Why Txn == Transcription?

Any other abbreviation is likely to work for "translation" as well and the x makes it sound a bit like transcription (txn would be pronounced something like "tixon").

I'm guessing tx came into existence first, and tx is closer to "transcripts".


text:
Java Protein Comparison Visualisation
I have project where I am required to build a visual analytics tool into an existing Java project to visually compare and analyse protein sequences from SwissProt, GenBank, EMBL etc. (downloaded from the databases in plain text and locally converted to XML)
Assuming I have the XML files of the protein sequences already, I am unsure regarding the following:

The best ways to visualize protein sequence differences.
(e.g. what types of graphs/ diagrams? Any sample images or links I could reference from?)

What information should be derivable from these visualizations and how they help researchers.

Any advice on the above?
Note: I have no background in bioinformatics.
Any help would be greatly appreciated!

jalview is one possibilty : www.jalview.org


text:
Getting Unique Identifier List for GEO Datasets NCBI
AIM: Download "Unique Identifier List" for the following query from GEO DataSets.
Query: ("Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism])) AND ("2020/01/01"[PDAT] : "3000"[PDAT])
which means, all RNASeq studies deposited on GEO in the year 2020 for humans, mice or rat.
Problem: I need the GSE ID list for ~9k datasets, but while trying to download the list of ids, it loads to a blank page and nothing happens. Also, clicking on "Next Page" gives error.
I have been trying for the last 3-4 days but it doesn't work.
Steps to generate file:
"Send To" -> "File" -> "Format" (Unique Identifier List) -> "Sort By" (Default Order) -> "Create File"

You can use Entrez Direct for this. The following returns Unique Identifiers which are just bare integers.
$ geo_query='"Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism]) AND ("2020/01/01"[PDAT] : "3000"[PDAT])'
$ esearch -db gds -query "$geo_query" | efetch -format uid > gds_results.txt
$ wc -l gds_results.txt
9981 gds_results.txt
$ head -n2 gds_results.txt
200134092
200120931

Instead, if you are looking for a way to get the GSE accessions, you can use the built-in xtract command to parse the XML returned by esummary as follows:
$ esearch -db gds -query "$geo_query" | esummary | xtract -pattern DocumentSummary -first Accession > gse_accs.txt > gse_accs.txt
$ wc -l gse_accs.txt
9981 gse_accs.txt
$ head -n2 gse_accs.txt
GSE165829
GSE165824

I think I would just try to do this with GEOquery.
https://bioconductor.org/packages/release/bioc/html/GEOquery.html

GEOquery would download datasets for a given accession ID. R provides a library https://www.bioconductor.org/packages/release/bioc/html/GEOmetadb.html which contains metadata for microarray studies. I couldn't find anything similar to this for RNASeq studies so had to resort to the GUI.

Oh you're right, I did something similar before, i though it was GEOquery but i guess i don't remember... Ok, can you download the text file like:
"Send To" -> "File" -> "Format" Summary (text)
And then do something like:
cat gds_result.txt | grep 'Accession: GSE' | cut -f 3 | cut -f 2 -d ' ' | sort > gse.txt


text:
BioPython - Retrieve sequence records from pubmed database
I've been assigned a task to fetch sequence records from "pubmed" database:

fetch 30 record of type "fasta" from (pubmed database, with term hemoglobin AND alpha) as a sequence records

It doesn't make any sense to me!
I know the PubMed database is for publications and citations, but how can I retrieve those records from it? is that even possible?
I tried to search for the IDs with a given term in it and then retrieve the records from another database (like "nuccore"), but I don't know if this is right.
Any ideas?

The PubMed database will not recover any sequence data - that's a database of publications. I think the question was Genbank

Answer from @devon-ryan, converted from comment:
They meant the nucleotide database.


text:
`TxDb.Hsapiens.UCSC.hg19.knownGene` and `EnsDb.Hsapiens.v75` can't be installed
I think there is a problem with TxDb.Hsapiens.UCSC.hg19.knownGene and  EnsDb.Hsapiens.v75.  I'm struggling to install all of these packages as they pause at  byte-compile and prepare package for lazy loading . Can anyone try to re-install them in their local machine because am not the only one who is having the problem, please? They used to work fine until a week ago!
As usual, I install the packages by:
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

I also tried installing them individually:
BiocManager::install("EnsDb.Hsapiens.v75")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

It pause at:
* installing *source* package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ ...
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading

sessionInfo( )
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] BiocManager_1.30.10 compiler_4.0.3      tools_4.0.3

I've tried with and without RStudio.
I'd be grateful if someone can try to replicate the error, please?

Can you edit your question to include how you are trying to install the packages? thanks.

Sorry @user438383, I forgot to add how I installed the pacakges. I edited the question now and added more details.

Thanks. Just tried it on macOs 10.15 and it seems to have installed OK. Does the same problem occur if you try and install a different package via bioconductor, or a package using install.packages()?

Oh, thank you. This is extremely strange because I tried it on my friend's computer and had the same problem. I can install any other package via Bioconductor without problem except for these three. Can I ask if you installed all three of them?

It's been a week and I've tried uninstalling R, RStudio and RFrameworks  countless times but still not working. I'm sorry but can I ask if if you used BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene")) to install them?

This sounds like a system issue of some sort possibly having to do with Big Sur. Have you tried to install any of these with conda?

I would vote to close this as an identical (towards content) question has been asked at Bioc support where this belongs https://support.bioconductor.org/p/p133849/ and there were also similar questions at Bioc towards these database packages just this week. it is New Year's eve so please give Bioc people time rather than opening questions at other communities, I guess this is something that they can and should answer.

Sorry, I only wanted to see if someone can replicate the error because am having it in two machines. My apologies again and feel free to close the question.

It just paused indefinitely?

Did you mean it also paused indefinitely in your local machine when you installed it? Thank you so much.

I am able to successfully install these packages using:
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

This completes without any errors / warnings (R v4.0.3 / macOS 10.16), so I don't believe the packages themselves/Bioconductor are 'broken', but I suspect you are having compiler issues due to the Big Sur update.
These are the steps I used to solve the problem:

Reinstall xcode command line tools (even if it says "up to date")

sudo rm -rf /Library/Developer/CommandLineTools
sudo xcode-select --install

Install gcc & llvm via Homebrew (instructions for installing Homebrew) or, if you already have gcc/llvm installed via Homebrew, skip to the next step

# This can take several hours
brew install gcc
brew install llvm

Once you have gcc & llvm installed via Homebrew:

brew cleanup
brew update
brew upgrade
brew reinstall gcc
brew reinstall llvm

Link some headers into /usr/local/include

sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/

# You can safely ignore warnings like this:
#ln: /usr/local/include//tcl.h: File exists
#ln: /usr/local/include//tclDecls.h: File exists
#ln: /usr/local/include//tclPlatDecls.h: File exists
#ln: /usr/local/include//tclTomMath.h: File exists
#ln: /usr/local/include//tclTomMathDecls.h: File exists
#ln: /usr/local/include//tk.h: File exists
#ln: /usr/local/include//tkDecls.h: File exists
#ln: /usr/local/include//tkPlatDecls.h: File exists

Create a new ~/.R/Makevars file (if you already have a ~/.R/Makevars file, save it in a different directory (away from ~/.R/)) and include only these lines in the file:

FLIBS=-L/usr/local/gfortran/lib/gcc/x86_64-apple-darwin19/10.2.0 -L/usr/local/gfortran/lib -lgfortran -lquadmath -lm
CXX1X=/usr/local/gfortran/bin/g++
CXX98=/usr/local/gfortran/bin/g++
CXX11=/usr/local/gfortran/bin/g++
CXX14=/usr/local/gfortran/bin/g++
CXX17=/usr/local/gfortran/bin/g++

LLVM_LOC = /usr/local/opt/llvm
CC=/usr/local/gfortran/bin/gcc -fopenmp
CXX=/usr/local/gfortran/bin/g++ -fopenmp
CFLAGS=-g -O3 -Wall -pedantic -std=gnu99 -mtune=native -pipe
CXXFLAGS=-g -O3 -Wall -pedantic -std=c++11 -mtune=native -pipe
LDFLAGS=-L/usr/local/opt/gettext/lib -L$(LLVM_LOC)/lib -Wl,-rpath,$(LLVM_LOC)/lib
CPPFLAGS=-I/usr/local/opt/gettext/include -I$(LLVM_LOC)/include -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include

Compile a package from source in R/Rstudio

# To check whether openmp is enabled, compile data.table:
install.packages("data.table", type = "source")

If your package fails to compile, a couple of SO users have had to install a fresh gfortran (re: https://stackoverflow.com/a/65334247/12957340), which you can download from https://github.com/fxcoudert/gfortran-for-macOS/releases/tag/10.2-bigsur-intel

Once you can install packages from source without error/issue, update Bioconductor
BiocManager::install(version = "3.12")

Then, finally, you should be able to install your packages using
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

You are amazing Jared, I'll go through the steps right now to see if it works, fingers crossed.

Awesome, it is working now. Thank you so much, Jared :)


text:
biopython-motifs , How can I create motif from sequences of different length?
I want to create motif of different length kmers.
Following code works for kmers of same length but not different length.What to do ?
from Bio import motifs

sequences = SeqIO.parse("myfile.fasta","fasta")
    instances = []
    for data in sequences:
      instance = data.seq
      instances.append(instance)

motif = motifs.create(instances)

Thank you for the edit @finswimmer

Motifs are not designed to be variable length in general since PSSMs/PVMs don't support insertions/deletions (indels). Instead, I suggest profile HMMs, such as those offered by HMMER3, which do support indels. You can parse HMMER3 output with Biopython via Bio.SearchIO.HmmerIO

All sequences represented in a Biopython motif object must be the same length. For example, a motif has a __len__ property, which has to consist of a single value.
There are other methods that do support variable-length motifs, I suggest that you look into them, but unfortunately biopython seems to be limited to fixed-length motifs.


text:
What are the applications of DNA or RNA pattern matching?
I'm assuming we don't do pattern matching in DNA or RNA for the fun of it. So I'd like to know what are the applications of pattern matching or where does it fit in in larger applications?
I'm a computer scientist and when studying certain topics in CS we see a lot of references to applying algorithms to areas of biology or genetics one of which is applying text pattern matching algorithms to DNA to find sub sequences. Whenever I see those kinds of examples I think "Oh that's pretty cool" but then I don't see how it's applied in the real world for benefit and I can't find anything to suggest where these algorithms are applied in these fields.

One useful application of pattern matching is for DNA sequence motifs. As an example of a direct application, there is a class of proteins called transcription factors (TFs), which bind to DNA and regulate the expression of genes. There are many TFs. Each TF likes to bind to its particular conserved sequence of DNA, with some wiggle room for some parts of the sequence.
Algorithms for finding and scoring these patterns (and mismatches in patterns) is important for identifying regions where regulation occurs, which proteins are involved, and for answering biological questions related to gene regulation and basic biological functionality — and malfunction, in case of disease.

These algorithms are used quite a lot, and I will touch only a few aspects:

Genetic sequence alignment. Different (but still closely related) organisms have very similar genetic sequences, which differ only at some positions due to point mutations, which can be nucleotide substitutions (another letter in a sequence) or insertion/deletion of a few nucleotides. Many biological questions such as studying new coronavirus strains, predicting genetic diseases, and identifying cancerous cells are based on such sequence comparisons. This requires aligning many sequences, using the appropriate algorithms. This may be not hard for two sequences, but for hundreds or thousands of them it may require serious computational effort.
Sequencing reads mapping Much of the current progress in genetics is due to the advent of the next generation sequencing, also called massively parallel sequencing, as opposed to Sanger sequencing used in the first human genome project (see here, for example). The basis of this technology is splitting the genetic sequence (a few billion nucleotides in a human) into small fragments of a few hundred nucleotides long and sequencing them simultaneously. One then faces the problem of putting them back together, which can be achieved either via matching them to a known sequence for a similar organism (sequence mapping) or by putting the reads together on the basis of their similarity (sequence assembly). There exist multiple sophisticated bioinformatic tools for achieving this (BWA, bowtie2, TopHat, Megahit to name just a few, so that you can google them).
Blasting Another application, closely related to my first point, is quickly identifying the sequence as belonging to a certain gene of a certain organism by comparing it to a large database. BLAST is one of the well-known collections of tools for doing this.
Just to mention (only a bit) more exotic applications: the algorithms associated with parsing context free grammars are used for RNA structure prediction.

Durbin's book is a good casual reading (for somebody with a computer science background) to get an overview without knowing too much biology.

The applications can be categorized based on whether you study the sequence for analysis, or for DNA synthesis and assembly.
In addition to TF patterns mentioned in this answer, other patterns for DNA binding proteins, enchancers, promoters, and the Kozak sequence (GCCGCC[A|G]CCATGG..) may also be useful for discovery of ORFs (genes) / regulation.
Certain methylase enzymes recognize a DNA pattern. For example EcoKI methylates adenine in AACNNNNNNGTGC (2nd A) and its reverse complement, GCACNNNNNNGTT (3rd A). Another such enzyme is EcoBI. This is useful for finding methylation sites, which may interfere with DNA restriction or expression.
Some restriction enzymes also recognize a DNA pattern: HinFI recognizes
5'GANTC
3'CTNAG

Other such enzymes are EcoP15I, EcoRII. Finding these sites are very important for designing a DNA cloning or manipulation experiment.
For DNA synthesis, it is useful to avoid homopolymers (N+) and repeated k-mers and hairpins (a sequence segment which has a reverse complement "nearby") as these can cause problems during synthesis.
Direct, tandem or inverted repeats can cause problems during cloning.
Finally, a perhaps trivial example is matching a polypeptide sequence, using DNA, as the genetic code is redundant. An example would be searching for target peptides in a protein.

To extend on the existing answers I would say that pattern matching algorithms are at the core of bioinformatics, certainly any NGS realted bioinformatics. The Smith-Waterman algorithm and its variants are probably the most used https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm .
Reference genome assembly (the human genome project for example) is based on kmer based sequence matching. Variant calling is based on mapping sequence to this reference, so thats pretty much all clinical bioinformatics. And a million other tools - one random example is dectecting horizontal gene transfer of antibiotic resistence genes.


text:
ATAC seq : Plotting insert sizes
My tutorial link
library(GenomicAlignments)

atacReads <- readGAlignmentPairs("MACS2_peak_call/ALL_ATAC_data/HSC1.bam", param = ScanBamParam(mapqFilter = 1,
            flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname",
        "mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))
# length(atacReads)
atacReads

atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)

library(magrittr)
library(dplyr)
library(ggplot2)
fragLenPlot <- table(insertSizes) %>% data.frame %>% dplyr::rename(InsertSize = insertSizes,
    Count = Freq) %>% mutate(InsertSize = as.numeric(as.vector(InsertSize)),
 Count = as.numeric(as.vector(Count))) %>% ggplot(aes(x = InsertSize, y = Count)) +
  geom_line()

fragLenPlot + theme_bw()

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

How to i specify the Grange to span the total chromosome. Now I use only "chr20" as shown in the example.

sorry don't quite get the question. is it about knowing the chromosome length from the bam file, or it's something to do with the plotting

Now Im too bit confused ,after reading paper where they don't mention which chromosome they taking and plotting.  It would be essentially this one " is it about knowing the chromosome length from the bam file" can i like  "chr1-chr20" and see the insert size distribution . I was not able to figure out how to define more that one chromosome in that Grange argument //

@StupidWolf for example in this paper https://www.nature.com/articles/s41597-019-0071-0#Sec1  Fig2(b).

You can exclude the which statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).
Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as bamPEFragmentSize from deepTools.

"bamPEFragmentSize" will see this

" you're seriously reinventing the wheel." yes since I didn;t read about your tool "bamPEFragmentSize" I used it , Works absolutely fine

is there any way by using this "bamPEFragmentSize" i can get smooth curve ?

Save the values to a file and replot however you like.

okay i will see the option to save the file ....right now i ran this "bamPEFragmentSize -hist fragmentSize.png --maxFragmentLength 100 -b HSC1.bam Mono1.bam LSC1.bam Blast1.bam --numberOfProcessors 10 --samplesLabel HSC Mono LSC Blast
"


text:
Gold standard benchmark
This page is claimed to contain a gold standard benchmark for viral genome assembly.
https://github.com/cbg-ethz/5-virus-mix
The claim is here:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/

Baaijens et al. De novo assembly of viral quasispecies using overlap graphs
Genome Research 2017 May; 27: 835–848.

"This data set was recently presented as a gold standard benchmark (Di Giallonardo et al. 2014)"
But at (Di Giallonardo et al. 2014) I can only find that mention:
"We previously estimated the frequencies of the virus strains in the 5-virus-mix by amplifying the protease gene using single-genome amplification (SGA), the current gold standard for studying diversity of virus populations".
That doesn't seem to mean precisely that the dataset is a gold standard benchmark for genome assembly but, as I understand, that SGA is a gold standard for something else.
My question is: is it? And who decided it is? Is still valid? Is there an official list of results somewhere for this benchmark set?
I tried to contact the owner of the page but still no answer.

Anyone can declare anything to be a "gold-standard", the term is meaningless. All if actually means is, "we think this is the test that everything else should be benchmarked against".

Thanks for raising the issue, it is certainly topical. The question is whether the data set they unit tested is reliable.
I think the issue they are trying get to is mosaic artefacts in the sequencing process and assembly thereon. This would be particularly true in a situation of de novo assembly which is what the authors are targeting. There are some slick methodologies for virus sequencing, which have been around for a while now, so in any case a 'gold-standard' would be outdated at present.
I do have to agree with you that the language particularly for the journal certainly could have phrased better. The data set was 2014 and genome circularisation (its a virus small genomes) to remove the otherwise high SNP error rate would have been around in 2017 when this was published.
I fully understand why the authors choose this language. A professional editor, not affiliated with the journal (without up to date technical insight), would likely recommend the same language.
I guess the question would be whether verification is achieved on other virus data sets. There is no shortage of them at present and its certainly an interesting tool.

Well, I was not asking that because of any suspicion about their results. My software (s-aligner) gets better results for this benchmark set and I just wanted to know how relevant is that. It also outperforms other software in many other benchmarks.

Thanks @juanjo75es this is a good subject to explore in a question. The algorithms are pivotal to understanding viral diversity. I put this as a separate question and notify you. I'll explain what I understand and why within the question.

@juanjo75es I have raised the matter in the following question, https://bioinformatics.stackexchange.com/questions/18747/metagenomics-pipeline-recommendations


text:
How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)?
Currently, I am implementing a descriptor of protein structure and would like to calculate properties per atoms like charges, hydrophobicity, hydrogen bond donors / acceptors, hydrophobicity, aromaticity etc.
I tried using PDB2PQR and succeeded with protein structures but it crashed while calculating nucleic acid bound proteins and choosing other force fields did not help.
Also, as far as I know, Rosetta tool can calculate hydrophobicity however, it assigns the property per residue but not per atom (correct me if I am wrong).
Any other suitable tool suggestions would be appreciated.

Per Atom
You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. OE1 and OE2 in glutamate behave identically).
So, the real question is whether you really care about the atom level in a classical mechanical setting?
Formal charges
A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say aspartate vs. aspartic acid.
PDB2PQR formats to a PQR file for the ligand docking application autodock. It does not assign how a residue is protonated.
In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.
Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.
A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (HIE in Amber) or delta nitrogen (HID), while it can be deprotonated (anion) or diprotonated (HID, cation). Rosetta switches between these two automatically.
Properties
You list a few properties of interest, which all seem to be based on the residue a priori of structure.
For example:

hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.
aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.
hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.

A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the "blur" of the atom. However, if you get the NZ atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse.
In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.
Here be dragons

Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.

However if for some exceptional reason this is required, say
figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:
import pandas as pd
import pyrosetta
pyrosetta.init()
pose = pyrosetta.pose_from_sequence("ELVISISALIVE")
target_resi = 9
residue = pose.residue(target_resi)
scores = []
for i in range(1, residue.natoms() + 1):
    iname = residue.atom_name(i)
    for r in range(1, pose.total_residue() + 1):
        other = pose.residue(r)
        for o in range(1, other.natoms() + 1):
            oname = other.atom_name(o)
            score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue,
                                                                                i,
                                                                                other,
                                                                                o,
                                                                                scorefxn)
        scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                         [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))

df = pd.DataFrame(scores)
df.pivot_table(index='atom_name_A',
               aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]


text:
What changed wwith Python in WDL change from draft 2 to version 1?
I'm updating my WDL from draft 2 (whichs works fine) to version 1. I'm stuck with a Py snippet
task getReadGroup {

    input {
        String inputFile
    }

    command {echo input = ${inputFile}}

    # command <<<
    # python <<CODE
    # print ("${inputFile}".strip().split('/'))
    # CODE
    # >>>

    output {
    Array[String] o = read_lines(stdout())
    }
}

This command gives the expected result:
input = /mnt/lustre/working/genomeinfo/data/20191015_LongRead/20191015_LongRead_1/data/LXBAB201271/LXBAB201271/LXBAB201271/20190417_0612_1-A7-D7_PAD44174_42dcde1a/fast5/test/

But the commented out python code:
print ("${inputFile}".strip().split('/'))

Gives an empty list - "${inputFile}" is empty in the code command for some reason!
Any ideas? I've tried using a File instead of String to make it just like the example in the docs (below) other than that I see no diff:
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

print (inputFile.strip().split('/')). ??

Is that a wild guess or do you know?

Don’t know what is WDL and how it works but it’s like the variable it’s inputFile , and $ goes with echo, but could be wrong. Can’t try by my own given I don’t know WDL

Given u don't know what WDL is it would be best if you didn't make suggestions...

Fine https://github.com/openwdl/wdl/blob/main/versions/1.0/SPEC.md#command-section got the use of placeholders    in whatever WDL is

What does print ("${inputFile}") alone gives you

Complements of Phil Shapiro;
Use ~{} instead of ${}.
So this example from the current docs is wrong:
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

It should be :
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("~{in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}


text:
Why did expression based subtypng of breast cancer gain much more acceptance than others
This is may not be entirely technical question but rather a academic question. But the technique behind the application is within the scope of bioinformatics. So I would try to ask here that:
In each cancer type, there have been tons of papers that describe RNA-Seq based subtyping resulting in molecularly distinct subgroups with survival difference. However, only very few of them gained acceptance at clinical setting, such as breast cancer (luminal subtype A, luminal subtype B, normal breast-like subtype, HER-2 overexpression subtype and basal-like subtype).
My question is: why did this molecular subtyping of breast cancer gain much acceptance but not others? Did this one own some advantages that other did not?

In two words: incidence and funding
I'm not an expert on this topic, but I assume it has something to do with the incidence of breast cancer itself:

Breast cancer is the most common cancer in American women, except for skin cancers. Currently, the average risk of a woman in the United States developing breast cancer sometime in her life is about 13%. This means there is a 1 in 8 chance she will develop breast cancer. This also means there is a 7 in 8 chance she will never have the disease.

This incidence most likely amplifies how big of a problem society/funders/governments consider breast cancer to be. For example, a study performed by Northwestern University’s Feinberg School of Medicine and Robert H. Lurie Comprehensive Cancer Center showed that breast cancer research is by and large one of the most well-funded areas of cancer research.

This higher amount of funding leads to more and better resources and research into every aspect of breast cancer, from identification to treatment. I would wager better funding has had the effect of increasing the level of research performed in the usage of sub-typing for identifying people at risk of getting breast cancer. Given positive results and lots of literature backing up such a method as useful, its levels of acceptance in clinical settings increased.


text:
Differential gene expression analysis of time series with replicates
I have a dataset that has two groups, perturbed vs control. Each group has 3 replicates. Each replicate has 8 time points. How do we do Differential gene expression analysis to find significantly expressed genes? Any packages in R or web tools etc?

R packages like DESeq2, edgeR, and sleuth will help you perform differential gene expression analysis.
These are typically used for control-vs-treatment experimental designs.
So you could use them like this and perform a bunch of pair-wise comparisons between the time points (i.e. what are the differentially expressed genes between timepoint $t_0$ and timepoint $t_1$, $t_0$ and $t_2$, etc.).
But all of these methods make use of generalized linear models to perform the statistical inference on coefficients.
For the simple control-vs-treatment design, above, this looks something like
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i X$$
where $g$ is the link function, $\mu_i$ is the mean expression of gene $i$, $\beta_i$ is the expression fold change estimate, and $X$ is the indicator variable saying if the sample is a treatment or control sample.
You can modify this equation by adding a time component to it.
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i^{(group)} X + \beta_i^{(time)} t$$
If you specify the timepoint of the sample in the experimental design matrix for DESeq2, edgeR, or Sleuth, you can get information about the effect of the perturbation as well as the effect over time.
You can try more complicated designs, too, like a coefficient for the interaction between timepoint and perturbation status, but with 3 replicates in each, it's probably best to stick with simple, linear models.


text:
Data format for pathway based clustering of samples
I came across this paper as one of the examples from this paper, this one
Figure 2. Host Protein Alterations in Infected iAT2s
(B) Dysregulated (FDR < 0.05) functional modules, including (bottom) U7 snRNA and HDAC-associated complexes.

So how should be the data frame should be for my clustering and visualization.
If i have to replicate the same for lets say 4 different time point how to do that.What i have done so far.
I have 4 condition C1,C2,C3,C4. I did pairwise comparison between C1vC2, C2vC3 , C3vC4.
So now I did pathway enrichment with the comparison above. I have for C1vC2, then C2vC3 and C3vC4.
Now I how to proceed what they have shown here. Do i intersect all the pathway that is common in all three? such that based on the any metric that is present I can cluster and show the differences or is there any simpler way.


text:
Load balancing numerous single-threaded processes with Snakemake on the grid
I'm familiar with running Snakemake in a DRMAA-enabled grid environment. (Most often, I write Python programs as wrappers and invoke the workflow via Snakemake's Python API.) Typically, I'm running workflows with rules that perform large, multi-threaded tasks. The rules are usually designed to utilize all available processors on a node, and only a single task is executed on any given node at a time.
But what if I have thousands of single-threaded tasks that I want to distribute among the requested grid nodes? If the corresponding rule(s) are marked threads: 1, is Snakemake smart enough to run multiple tasks on each node? How does it know how many tasks can be run on an individual node—or how does the user specify a limit?


text:
How to represent random sequence elements in SBOL?
In our lab, we work with synthetic biology components with partially random sequences (similar to work in directed evolution). So, for example, we have a plasmid design with several components, including one protein coding sequence that will be created from a wild-type sequence via mutagenic PCR.
At the design stage, we don't know the exact DNA sequence for any of the resulting variants of that protein coding sequence. We just know the wild-type sequence and the average number of mutations we'd like to apply.
How can we represent that design in SBOL (via the Component class)?
Then, when we implement the design, we end up with a single sample that contains plasmids with 10^5-10^6 (or more) different protein coding sequences. We then measure the complete DNA sequence for every plasmid variant, but I don't want to end up with >10^5 different SBOL definitions to describe the sample.
How can we best represent that sample in SBOL (via the Implementation class)?
Finally, we also end up measuring the performance of every plasmid variant, and we identify the specific protein coding sequences for useful and/or interesting variants. We then design and build specific, clonal plasmids with those useful/interesting protein coding sequences.
Representing those clonal plasmids as SBOL Components seems pretty straight-forward, but is there a good way to capture the relationships with the original random-library plasmid design and the specific implementation of that design that we built and measured?

This sounds like a job for the CombinatorialDerivation class, which is designed for representing libraries and the relationship between derived designs and the originals.
About CombinatorialDerivation
In SBOL 3, a CombinatorialDerivation takes a Component in its template property to represent the system that you intend to vary:

Some of the Features of the template Component are designated as the variables of the library via VariableFeature child objects of the CombinatorialDerivation, which point at the variables with their variable property.
The values that can be filled in for each variable are then designated via the variant* properties of the VariableFeature.
Finally, the distribution of values is determined by the strategy property of the CombinatorialDerivation and the cardinality properties of each variable.

Representing Random Mutagenesis
In the case of random mutagenesis, if I'm understanding your design correctly, you would represent it with a CombinatorialDerivation as follows:

The template is Component of type DNA with the original sequence. Each site designated for mutagenesis can be designated with a SequenceFeature.
For each of the mutagenesis sites, create a VariableFeature in the CombinatorialDerivation with the corresponding SequenceFeature as its variable and a set of variant values for each of the mutant options set by your primers. The cardinality is probably one (if every site has to be hit by something in your protocol) or zeroOrOne (if sites might or might not be hit).
The strategy for the CombinatorialDerivation is sample, since only a few of the many possibilities will actually be realized.

Representing Samples and Characterized Constructs
Now, when you take it to the lab, you can represent your aliquot with an Implementation whose wasDerivedFrom points to the CombinatorialDerivation, representing the fact that you've got a randomly generated library whose contents are not precisely known.
Finally, once you've got your isolated products represented as Components, there is a recommended best practice to use wasDerivedFrom links to connect the derived Component and its varied Features back to the CombinatorialDerivation and its template. If you make a Collection of your results, the Collection should also link to the CombinatorialDerivation.
This should then capture the relationship between all of the products, their associated data, and the original mutagenesis design.

I have some code for doing a similar task that I'm happy to share if that's of interest; it happens to not have yet been posted publicly yet, but I can quickly clean it up a bit and do so if desired.

No need to rush posting of the code. I'm working toward the framework outlined in my question, but have a lot of other things to do first. I just wanted to make sure I was pointed in the right direction.

I do have one clarifying question/comment: From the "Representing Random Mutagenesis" part of your answer, if we are going to use mutagenic PCR to make the protein variants, then every base in the protein coding sequence would be a potential mutagenesis site. So, I would need to define ~1000 SequenceFeatures, one for each base in the protein coding sequence. Is that right?
This seems a bit cumbersome, but but not too hard to write code for. Does the idea of a SequenceFeature for each base in a protein coding sequence raise any red flags for you though?

@DavidRoss I think I misunderstood your protocol: I thought it was something like this protocol for targeted swaps with a random subset of a fixed set of variations.  If you're randomly varying every base in a large region, then making a SequenceFeature for every base in the sequence does raise a red flag to me.  I don't think we've worked out a best practice for this case yet, but my thought would be to mark the whole region as one SequenceFeature and tag it with a variant of all NNNNN.... bases

That way, you're saying "replace this region with indeterminate bases" and any variant you find will be an appropriate realization of that specification.


text:
How to obtain NCBI Nucleotide (Entrez) IDs for all E coli sequences?
I would like to gather a list of ALL E coli NCBI (Entrez/Nucleotide) IDs so that I can download fna and faa files, through a function like this:
# python code #
from Bio import Entrez
from glob import glob

def dl_genome(id, folder='genomes'):
    files=glob('%s/*.fna'%folder)
    out_file = '%s/%s.fna'%(folder, id)

    if out_file in files:
        print (out_file, 'already downloaded')
        return
    else:
        print ('downloading %s from NCBI'%id)

    Entrez.email = ""
    handle = Entrez.efetch(db="nucleotide", id=id, rettype="fna", retmode="text")
    with open(out_file, 'w') as fout:
        fout.write(handle.read())

for strain in targetStrainIDs:
    dl_genome(strain, folder='genomes')

My problem is that I don't have a good way to get the E coli strain names and their corresponding NCBI IDs for Entrez. Is there a way to get these NCBI IDs without manually searching for them?
Context
I have been trying to develop an E coli pan-genome for one of my projects. I initially just downloaded the relevant genomes (fna, faa, & features.tab files) from the PATRIC database. However, I've come to realize that the reference sequence on PATRIC has issues. Specifically, the reference sequence on PATRIC's fna/faa files for E coli K12 MG1655 notably do NOT list genes like b2031 or b2033, despite both being present in this strain.
The NCBI reference genome DOES include these genes (with proper annotations), so I'd like to just download what I need from NCBI instead. However, I'd like a systematic way to download all the fna/faa files since my pan-genome is eventually going to be built from hundreds of of these fna/faa files.
Any help would be appreciated.

You should check out NCBI Datasets for this sort of thing. Since this is bacteria, you won't be able to use the web interface but there's an API and a command-line tool that you can use.
For example, you can use the datasets command line tool to download the E. coli RefSeq genome assemblies as follows:
## see how many genomes are in scope
$ datasets summary genome taxon 562 --limit none --refseq
{"total_count":22287}

## download data excluding gff3, rna and genome seq files
## as an example but drop --exclude-seq to download genome seq as well
$ datasets download genome taxon 562 --refseq --exclude-gff3 --exclude-rna --exclude-seq --filename ecoli_genomes.zip

Do you mean the NCBI taxonomy IDs? If so, I've personally always preferred the way UniProt displays NCBI taxonomy, and it's easier to download IDs in mass.
For example, if you find E.coli's entry: https://www.uniprot.org/taxonomy/562
You can then click on "all lower taxonomy nodes", which will show all E.coli strains. If you click 'Preview first 10' before downloading, you'll see that the first column is the NCBI taxonomy ID of each strain of E.coli.
For example, E.coli strain 0.1197 is reported to have the ID 1005443 from this preview page. Indeed, if we search this ID on NCBI's taxonomy page, we do get the correct entry: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1005443
If this is what you want, click 'download all' to get a nice tab-separated file of all the strains of E.coli and their NCBI taxonomy IDs. This file could then be read and processed before you send API requests to Entrez.
Hope this helps!


text:
How to get 2d/3d structural coordinates from a PDB file?
Is there a way we can get 2d/3d structural coordinates like in (x/y/z) manner from a pdb file.
I have looked a believe that the coordinates are in below line:
ATOM   1067  NH1 ARG A 141     -10.147   7.455  -6.079  1.00 23.24           N

but i am unable to make a sense from it as i don't know which one is (x/y/z) in (2d/3d) coordinate system.
If anyone describe it, that'll be a great help.
Thanks

Molecular structures are 3D.
The records (lines) in the PDB files are in so-called fixed column format.
The fields don't need to be separated with blank characters, but they must be in specified columns. For example, the X coordinate is in columns 31-38.
Here is the specification of the ATOM record from PDB v3.3.
COLUMNS        DATA  TYPE    FIELD        DEFINITION
-------------------------------------------------------------------------------------
 1 -  6        Record name   "ATOM  "
 7 - 11        Integer       serial       Atom  serial number.
13 - 16        Atom          name         Atom name.
17             Character     altLoc       Alternate location indicator.
18 - 20        Residue name  resName      Residue name.
22             Character     chainID      Chain identifier.
23 - 26        Integer       resSeq       Residue sequence number.
27             AChar         iCode        Code for insertion of residues.
31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
55 - 60        Real(6.2)     occupancy    Occupancy.
61 - 66        Real(6.2)     tempFactor   Temperature  factor.
77 - 78        LString(2)    element      Element symbol, right-justified.
79 - 80        LString(2)    charge       Charge  on the atom.


text:
Input from multiple channel, merge and use in nextflow
I need to merge multiple channel to reuse them in two processes. So far I tried as follows:
channel.from([15,20,21]).into {k1;k2}
channel.from(['a','b','c']).into {m1;m2}
channel.from([1,2,3]).into {l1;l2}
para1=m1.merge(k1).merge(l1)
para2=m2.merge(k2).merge(l2)
process p1{
input:
each m, k, l from para1
"""
echo $m, $k, $l
"""
}
process p2{
input:
each m, k, l from para2
"""
echo $m, $k, $l
"""
}

The error I am getting:
Unknown process directive: '_in_each'

Did you mean of these?
    each

If I use para1.view() I can see the merged elements as expected. Any help

I don't think you can unpack tuples using the each qualifier like that. You can, however, unpack them in the script block:
process test {

    input:
    each triple from para1

    script:
    def (m, k, l) = triple

    """
    echo "${m} ${k} ${l}"
    """
}

Consider also that the merge operator is deprecated and will be removed from future Nextflow releases. Instead, you may need to join by index, for example:
def indexedChannel( items ) {
    return Channel.from( items.withIndex() ).map { item, idx -> tuple( idx, item ) }
}

k_ch = indexedChannel( [ 15, 20, 21 ] )
m_ch = indexedChannel( [ 'a', 'b', 'c' ] )
l_ch = indexedChannel( [ 1, 2, 3 ] )

m_ch
    .join( k_ch )
    .join( l_ch )
    .map { idx, m, k, l -> tuple( m, k, l ) }
    .into { para1; para2 }

Thanks @Steve. Personally I liked merge option (fewer lines). The processes ran without any error but how could I see the outputs?

No worries @zillurrahman. The test process above will just echo to stdout and this should be captured in the '.command.out' files in the process work dir. HTH.


text:
How to draw a 3d structure using 3d coordinates from pdb file?
I am trying to create a 3d structure from 3d coordinated from a pdb file.
I am trying to draw 3d coordinates from pdb file "ATOM 1 N VAL A 14 22.692 -27.591 -39.402" and using webgl javascript.
I've managed to get the 3d coordinates from a pdb file but i don't know the exact logic that is used to create 3d structure from the coordinates found in pdb file.
I tried those coordinates and managed to create a cloud of 3d points but it is not anywhere near the structure that is found on pdb site.
So wondering how the structure is created on the pdb website?

can you please add the code you tried so far?

JSmol sounds like what you want - it'll let you display 3D structures from PDB files onto webpages using JavaScript. If you just want to view the structure rather than display them on webpages, there's many programs that exist already like Swiss-PdbViewer and Pymol. There's even web-interfaces that let you upload a PDB file and they display it for you like Jena3D Viewer.

If you want to share a PDB file, you could try Michelanglo. It allows you to upload a PDB file (among other things) edit a description panel (which can feature special links that control the protein view and representation) and share the link with whomever —without needing any installations on either side.
Here is the page of shared protein: gallery
(disclaimer: I am the admin).
If you are interested in implementing a molecular viewer on your site, I'd suggest NGL. The PDB website uses NGL too, but is recently switching to MOL* ("molstar"), which is a new development, so less documented, but fixes some corner cases and is marginally faster and is a collaboration with a PDBe. NGL does not depend on any technology that will be depracated mind you, so it actually is a very solid choice for a JS molecular viewer.
The PDB format is a working format, not a deposition format, as is mmCIF. As a result NGL and Mol* have some issues when it comes to PDB files, such as lack of secondary structure records (HELIX/SHEET) results in anaemic sheets, fancier linkages (LINK/SSBOND reconds) are not handled only CONECT records.

disclaimer: I am the admin ..... I though you were the creator (!!) of Michelangelo ?

@pippo1980 Ops. Yes, I am: I was trying to be subtle, but at the same time abiding by the transparency rules/guidelines/custom of SO.

I had a look at Michelanglo and as a third  party I can say that it rocks !!!! Very nice tool indeed

Thank you —that is always nice to hear!

As others said, there already exist javascript libraries to visualize a protein. The most complete and recent one would probably be Mol*, which is what is used on RCSB PDB website. Mol* actually works using NGL, a WebGL engine. But probably the oldest library used on RCSB PDB to visualize proteins would be JSMol. You can actually switch between these viewers on RCSB PDB 3D view using "Select a different viewer" in the bottom right (try it out on 1UBQ)
Now to address the first exercise: there are several things going on after getting the atom coordinates. A PDB file may or may not contain a CONECT part, where bonds between atoms will be clearly stated (contrarily to a CIF file which I believe always contain the bonds). Henceforth these packages may or may not have to guess atom bonds. Given that a PDB file also contains the residues of each atom it's fairly easy to do, less so when dealing with small molecules.
However what you see most of the time is a cartoon representation of proteins, where each secondary structure is represented differently. The library will compute the secondary structure of each residue based on the dihedral angles of the backbone, then represent them. As of why they look like this, check the work of Irving Geis and Jane Richardson.
So if you want to have a website showing protein structures I would strongly advise to use an already existing library, although the exercise of coding your own can always be interesting. If you want to visualize structures on your own computer without using a web interface I can only recommend PyMOL, VMD and ChimeraX.

I am not good as you are and needed to use a parser for the pdb file
see https://mmcif.wwpdb.org/docs/software-resources.html:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  7 17:40:03 2020

@author: Pietro

"""

import matplotlib.pyplot as plt  #matplotlib

from pdbx.reader.PdbxReader import PdbxReader #PDB file parser

import numpy as np #numpy

openz = open('./1msc.cif')

print('opened')

# Initialize a PdbxReader object with the input file handle

pRd = PdbxReader(openz)

# Initialize a list to be propagated with DataContainer (and/or DefinitionContainer) objects (of the DataContainer class, which inherits from ContainerBase) parsed from the CIF file, where data blocks map to DataContainer objects

data = []

print(data)

# Call the read(self, containerList) method with your list

pRd.read(data)

print ('data ###################################################################\n' , data)

print('data lenght' , len(data))

# Your list is now propagated with one or more DataContainer objects, which represent data blocks. To get the first data block, just use list notation:

block = data[0]

print(' block ##############################################################\n',block)

pippo=block.getObjNameList()

print(pippo)

# To retrieve a category object, use the getObj(self, name) method

atomsite = block.getObj("atom_site")

print('atom site : ', atomsite)

atom = atomsite.getValue("group_PDB",0)

# print('atom sites id :', atom)

atomlist=[]

i=0
while True:
    atom =  atomsite.getValue("group_PDB",i)
    atomid = atomsite.getValue('label_atom_id',i)
    if atom =='ATOM':
        # print(atomsite.getValue("group_PDB",i))
        if atomid == 'CA':                         #select only carbon alpha (CA) to work with

            atomlist.append([float(atomsite.getValue('Cartn_x',i)),\
                           float(atomsite.getValue('Cartn_y',i)),float(atomsite.getValue('Cartn_z',i))])

        i +=1
    else:
        print('loop ends')
        False
        break

atomlist=np.array(atomlist)

print('end')

print(atomlist)

print('lunghezza : ', len(atomlist))

fig = plt.figure()
ax = plt.axes(projection="3d")

for p in atomlist:       #print CA as dots
    ax.scatter3D(p[0], p[1], p[2])
    print(p)

for i in range(0,len(atomlist)-1):   #connect with a line CA n with n+1
    x,y,z = [atomlist[i+1][0],atomlist[i][0]],[atomlist[i+1][1],atomlist[i][1]],[atomlist[i+1][2],atomlist[i][2]]
    ax.plot(x,y,z, color='red')

plt.show()

my input is https://files.rcsb.org/download/1MSC.cif
Result is just the CAs (carbon alpha) shown as dots and connected to form the main chain:

Didn't manage how to show the dimer though