text:
Duplicated genes and genome fragmentation
This is a follow up question to this post: How can I find out if my gene of interest is duplicated?
Briefly, to identify homologous genes, I performed tBLASTn analyses against my genome of interest using various sequences of a closely related species as queries for the BLAST analyses. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome. These genes are involved in different organelle functions (i.e., they are not surface antigens) and I want to find out whether these sequences represent duplication events or not. These analyses are being performed on protozoans parasites that belong to the same genus. The query genome is properly annotated, but the target genome is not (it is a whole genome shotgun project that has been assembled at the scaffold level).
Following the advice of @Michael G., I re-analyzed my sequences and blasted them against the entirety of Genbank. If there are two sequences for the same protein (for e.g. protein 3-1 and protein 3-2) and I compare them among themselves they are very similar (over 70% identity at the protein level). When I blast each sequence against the whole protein RefSeq database, the first hit is my query and the rest of the hits correspond to the same protein in other organisms. I believe that these sequences represent duplication events and it is possible that they are paralogues, but more evidence is needed to classify them as such.
However, when there are more than two sequences for the same protein, I found something different. For e.g., say I have three sequences for the same protein: sequence 5-1, sequence 5-2 and sequence 5-3. Sequences 5-2 and 5-3 are very similar (78% identity at the protein level) but when I compare sequences 5-1 and 5-2 I find that they are overlap in the middle, i.e., they are part of the same sequence. Thus, it appears that I don't really have three sequences for this protein but two. Blastp analyses against the whole RefSeq database give the expected result (first hit is my query sequence, the rest of the hits are same sequence in other organisms). So, once again I have what appears to be a duplication event, but in this case one of the sequences appears to be split in two -- could this be because the target genome is not completely assembled? Why could this be happening?
Thank you very much for your help.
-Leah
Edit: The new results are posted in this thread: Duplicated genes and genome fragmentation - Part 2

Thanks @leah, let me think. The "duplication event" will depend on the sample of other protist genomes available.

Thanks for your reply @MichaelG. I am sorry, but I don't quite understand... What do you mean by "the sample of other protist genomes available"?

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this? There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

Can you two please convert this to a Q&A instead of a discussion thread? Leah, please take all the information and make it into one, comprehensive question. @MichaelG. please post one answer that answers that comprehensive question. I don't really know what else to do here without deleting everything. Please don't use this site like a forum.

In the mean time, I have closed this as too broad since it is very much too broad for this site in its current form.

Hi @terdon, the issue is that it is different questions each time. The current question regards long branches and long branch attraction (LBA) in a phylogeny. This could be posted as a separate question. What is a Q&A thread?

@MichaelG. the idea here is to have one question and then everything else would be an answer to that question. If the question requires discussion like this, then it is too broad for the SE model and should be closed, and the question posted on SeqAnswers or BioStars instead since those to discussion. SE tries to have a high signal to noise ratio, which means that questions and answers need to be clearly separate and easy to find. When you have multiple back and forths like this, everything is confusing.

[cont] For instance, because of the voting system, there is no guarantee that the order of the answers will remain the same, so it's impossible to figure out who's answering whom and which answer is posted in answer to which question. The whole thing becomes a mess and ends up helping nobody except perhaps the people involved in the discussion. But we want to set up a resource that will help everyone, including the user who stumbles upon this question in 5 years time where the order of the various answers will be totally different. :/

In summary, apologies for the rolling nature of the question and answer and I will know in future. It is an interesting problem. If the OP targets specific questions, thus 1) long branch attraction (which is a big topic), 2) state-of-the-art tree amino acid tree construction (which is applicable to other studies), it should be compliant with forum rules?

The strategy from here is to determine duplication/parology through phylogeny, notably amino acid phylogeny using statistical support.
The basic result is whether the genes under investigation are:

A. monophyletic to the given species under investigation
B. Paraphyletic against protist species which are not under investigation,

Point A is compatible with paraology/duplication within the species. Point B infers the duplication is ancient.
The skill/difficulty is building a comprehensive alignment.
Before we start,

Parology is usually underestimated within a genome because repetitive sequences are often collapsed in the assembly process. Thus a de novo assembly around the gene of interest in both annotated and shot-gun assemblies might yield better results, ie. more duplicates. You have already seen what appear to be N-terminal and C-terminal ends of a related gene fragmented in the assembly process.
The tBlastn results looked quite cool and need to be included.

The core approach is to blastp your query sequence against the full NBCI database. Blast has an option to automatically download all hits (or all selected hits), e.g. in fasta format. This will provide a multi-species data set to understand the evolution of the protein of interest. This feature of blast is very cool and saves a lot of work. You want to combine the output of point 1 and point 2 in the fasta blastp list and align everything via ClustalO or Muscle (any alignment program). Then manually edit the AA alignment via an alignment editor (I use Seqotron [Mac]) and thereon produce a very basic tree using neighbor-joining just to check everything looks ok. The NCBI blast trick might need a better explanation so let me know if you have trouble figuring this.
If you post the tree as a PHYLOGRAM (I don't need to see the sequence names - just the shape and branch lengths) we can then see whether it looks reasonable before commencing a heavy duty algorhitm to identify support for the competeting hypotheses.
It is not ideal doing this via a message server, but its doable.

Thank you very much for the detailed explanation! I understand the basics of what you described, so I will give it a try. I hope that I will be able to download all the hits from blast... I will let you know.  Again, many thanks for all your help.

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?
There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

So this is what I did to produce the phylogenetic tree:

I performed a BLASTp using the query sequence (lets call it sequence 5 from species A, which is thoroughly annotated) against full NCBI protein database (ref-seq database)
I downloaded all the hits in fasta format (very cool feature indeed!)
Here I was a bit confused, so what I did was to take the three hits I got in my target genome (sequences 5-1, 5-2 and 5-3 from species B, which is poorly annotated) and added them to the fasta blastp list
I took the fasta blastp list and loaded it into Seqotron.
I aligned everything using the Compute - Alignment - MUSCLE option -- I didn't do any manual AA alignment (not too sure about this part)
I produced the tree using Compute - Phylogenetic Tree (default parameters: distance matrix: K83, resampling: none).
This is the tree I got. Sequences from the target genome (species B) are shown with a red rectangle:

Please let me know if I missed something or if I made a mistake somewhere. Many thanks!
-Leah

Cool @leah, this looks good. There is an issue with 2 of the taxa, the ones in red, otherwise if this was correctly rooted it looks fine. I'll post further things to check much later today

Many thanks for your response!

What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?

Yes, I have performed a blast back over the original highly annotated genome and, as expected, the first hit is the sequence that I used as query in the first blast analysis.

There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

I am working with Amoebozoans -- specifically, with members of the group Discosea - Centramoebida.
Once again, thank you very much for your help.

That is fine. From here it's difficult stuff and you need to build trees, providing advice over a forum like this is tricky but we'll have a go.

The first thing is the Sequence 5-2 and 5-3 are on very long branches. This is often indicative of sequence problems. You mentioned 5-2 and 5-3 are possibly two fragments of the same gene. My advice is to do a psuedo-annotation and reassess the phylogeny to see whether that super long branch length falls back in line
The branches close to 5-1 are of interest and their genomes should be further investigated.
You need to find the root of the phylogeny using external criteria. I can't do this because I've no idea was the species are.

Long branches are bad news in phylogeny, particularly if artefactual because of a phenomena known as 'long branch attraction'. Don't use your tree to understand duplication until this is resolved etc ... because the results can be distorted by LBA. Thus you must either delete them from the tree or perform a psuedo-annotation and redo the tree.
I was reluctant to forward the idea of phylogeny, because it is not just a bit of coding but an entire dialogue. Anyway by correcting the presumed misannotation you have sufficient results to move the investigation forward. You cant delineate your hypotheses until the root is correctly determined and then you set this as your outgroup under a more robust tree-building algorhithm.
I think we are out of space on this thread (you're allowed 1 post, I'm allowed 2), but clearly things have moved forward. Again, please note once the long-branch is gone you can end up with a different tree.

Thanks for the detailed explanation. I will annotate sequences 5-2 and 5-3 and run the analysis again to see if the branch length is shortened. I will also look into the branches close to sequence 5-1 and see what I can make of them. I will post those results in a new thread. I will also look into finding the root of the phylogeny, but I am not so sure how to go about it... Which external criteria could I use? Thanks again for everything!

The new results are posted in this thread: https://bioinformatics.stackexchange.com/questions/10365/duplicated-genes-and-genome-fragmentation-part-2


text:
How to run function "convolve" in parallel in R?
I have a vector of 200 million numbers. It has been running for over 24 hours and still not done.
x <- data   # the target sequence
k <- exp(-(-win:win)^2/(2*sig^2))/(sig*sqrt(2*pi))  # gaussian kernel
ptm <- proc.time()
y <- convolve(x,k,type="open")
proc.time()-ptm

How to parallelize the convolve function?
To parallelize, I tried splitting up the vector into pieces:
> x <- data[1:1e7]
> ptm <- proc.time()
> y <- convolve(x,k,type="open") # smoothed sequence
> proc.time()-ptm
   user  system elapsed
 40.743   0.206  40.940
> x <- data[1e7+1:2e7]
> ptm <- proc.time()
> y <- convolve(x,k,type="open") # smoothed sequence
> proc.time()-ptm
   user  system elapsed
320.193   0.708 320.871

However, splitting it up and then adding back the convolutions will be tricky.
Does anyone have other ideas?
This data is based on biological data. We want to cluster values in the vector based on nearby vector indices.
Reproducible example with 20,000,000 random numbers should finish in about 5.5-6 minutes. But if you try running 200,000,000 locally or on HPC, it will take a long time (my HPC has limit 99 hours and it timed out at 99 hours unfinished).
win <- 100        # window length of gaussian kernel (units acc to orig sequence)
sig <- 5        # std deviation of gaussian
k <- exp(-(-win:win)^2/(2*sig^2))/(sig*sqrt(2*pi))  # gaussian kernel
l <- 20000000
x <- round(runif(l, 0.0, 1.0),1)   # the target sequence with peaks
t <- seq(0,1-1/length(x),1/length(x))
ptm <- proc.time()
y <- convolve(x,k,type="open")[(win+1):(win+length(x))] # smoothed sequence
proc.time()-ptm

Is the function vectorized? Did you try the same function with less numbers? If you need to make a cluster, why do you need to make the convolution? Perhaps you will get the attention of those that could help you better on stackoverflow.com.

Hi IIrs: The convolve function is a function in R. The input data (x <- data) is a vector. We need to group based on nearby indices in the vector, so we cannot use other clustering methods that were built, for example, clustering gene expression. I can re-post in stackoverflow, but I have heard re-posting the same question is not good (and I am new to the stackoverflow community). I can try re-posting though. Thank you.

What kind of vector is it? A list or a numeric one? Could you please post a minimal reproducible example (with faked data if needed)?

Hi IIrs, It is a numeric vector. I've troubleshooted and number of decimal places does not affect time to finish. I have edited the post above with a reproducible example with 20,000,000 random numbers. The time to finish with random numbers matches the time to finish with subsetted real data.

I don't know what is supposed to do but maybe the problem is that there isn't any continuous periodicity to use fft on it??

source: Stack Exchange

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

sort: votes


text:
Create a Heatmap that uses modules of genes rather than individual genes in Seurat
Is it possible to create a Heatmap in Seurat that takes in and displays modules of genes as its features rather than individual genes?  So if I have a module of genes associated with a trait or phenotype I can compare the expression of that module across clusters against the expression of other modules of genes.
Currently the way that I am trying to do it is similar to the code below (In this example I am only using one module to test if it works):
Module_28 <- readLines(con = "Path_to_file/module28.txt")
Module_28 <- list(Module_28)
Combined <- AddModuleScore(Combined, features = Module_28, ctrl = 5, name = 'Module_28')
names(x = Combined[[]])
avg <- AverageExpression(Combined, features = NULL, add.ident = NULL, return.seurat = TRUE, verbose = TRUE)
DoHeatmap(avg, features = c("Module_281"), group.by = "ident", size = 3, angle = 0, draw.lines = FALSE) + theme(axis.text.x = element_text(size = 7, angle = 45), axis.line = element_line(color = "#ffffff")) + scale_fill_gradient2(low = '#1000ff', mid = "#ffffff", hight = '#aa0101', space = "Lab", na.value = "#ffffff", midpoint = 0, guide = "colourbar", aesthetics = "fill") + coord_flip(expand = TRUE, clip = "on")

Often When I try to run the DoHeatmap function I get the following error:
Error in DoHeatmap(avg, features = c("Module_281"), group.by = "ident",  :
  No requested features found in the scale.data slot for the RNA assay.

What's wrong with your current approach?

I get an error that I'll add to my post above so you can see it.

Did you scale your data at some point?

@haci yes I did scale the data, with the ScaleData() function.

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


text:
DMPfold: Understand the run of DMPfold "run_dmpfold.sh"?
I ran dmpfold using this command:
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

It produced files like this:
sequence_sub_embed_1.pdb ..... sequence_sub_embed_500.pdb
And, now I see some more files like this being produced:
sequence_1.pdb .... sequence_2.pdb .....etc
Can you please let me know what is dmpfold doing here? And can I stop dmpfold in between when desired number of models are generated? And, while it is running, can I move the previously made pdbs to new location becuase of disk space issues or all pdbs have to be in the "dmpfold-lysin" as stated in the command-line above? Thanks.

From the docs:

The final model is final_1.pdb and other structures may or may not be generated as final_2.pdb to final_5.pdb if they are significantly different. Many other files are generated totalling around 100 MB - these should be deleted to save disk space if you are running DMPfold on many sequences.

These other files include intermediate CNS files (sequence_sub_embed_1.pdb etc.), output files before clustering (sequence_1.pdb etc.) and saved ensembles of previous iterations (ensemble.1.pdb etc.).
Let the run finish or you won't get the final output files. If you want to use all the generated structures, use the ensemble.n.pdb file from the last iteration (i.e. the highest number). You can also use the other ensemble files if you want to include models from previous iterations.
I would strongly recommend using default "3 50" parameters until you get things working. Don't move files during the run. You really don't get much benefit of going beyond 50 models per iteration, trust me. We tried it relatively extensively!

Just one question: "output files before clustering" and "ensemble files": are these also 500 + 500 in number?

Yes, after a run has finished the sequence_n.pdb files will be 500 in number and correspond to the 500 structures in ensemble.n.pdb. Other ensemble files contain structures from previous iterations.

I noticed a behavior of the program I have: I left "run_dmpfold.sh" to run for two days without touching the system using the command line I mentioned, at that time dmpfold had produced  "sequence_71.pdb" file, and after 2 days when I check, there were no more number pdb files than "sequence_71.pdb". When I logged into the system as it was on sleep mode, the number of "sequence_71.pdb" stared increasing to "sequence_72.pdb"..etc. Hence, actually 2 days were lost where dmpfold was not producing files. Can you please let me know why this might have happened?

It sounds like your machine went to sleep, which will have paused all running processes. This isn't a DMPfold issue, it's more an OS issue. Have a look at tools like https://mac-how-to.gadgethacks.com/how-to/take-control-your-macs-sleep-functions-with-these-commands-0168109/ to sort this out.

yes, I corrected it by "aborting suspend" in my system settings. Thanks

Another query to clarify, I see that after generating "sequence_n.pdb" files, these files are no more in the folder? And, now I see "ensemble_1.pdb", will there be more ensemble files like "ensemble_500.pdb"? I mean when can I expect "final_1.pdb' file? Also, I see files like "CLUSTERS_001.pdb". thanks!

The sequence_n.pdb files are temporary and are combined into an ensemble.i.pdb file after each iteration i. That file has many structures in. After the run is complete there will be a final_1.pdb file. Don't worry about the cluster files. If you are interested in what they do, consult the run_dmpfold.sh script.


text:
Dmpfold error "No such file or directory: 'best_qdope.pdb'" when running "run_dmpfold.sh"?
I am running DMPfold, "seq2maps.sh" ran successfully, but whrn running "run_dmpfold.sh", it gives me this error:
   Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
seed =  30603
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/run_dmpfold.sh: line 177: c600m2-tail-before-lysin-sequence_[1-9]*.pdb: No such file or directory
No atoms read from PDB file
0 models read from PDB file
Outliers:

Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 82, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 55, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 75, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 53, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
No valid current contacts.

Can you please let me know how this can be solved?
the comand line is:
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

Thanks!

Please [edit] your question and show us line 177 of /media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/run_dmpfold.sh. As you can see in the error message, the script is looking for files named c600m2-tail-before-lysin-sequence_[1-9]*.pdb and it is doing so on line 177. What is this script? Where did you get it? Could it have hard-coded input files? All of your errors are about missing files.

The target variable, which is set to the first part of the sequence file name, should be sequence according to your given command line. However the error message thinks it is c600m2-tail-before-lysin-sequence. Have you edited the command line for simplicity?
I'm a little worried about the floating point errors, which are almost certainly from CNS as there are 3 CNS calls and the error shows up 3 times. I would test your CNS installation to see if everything is okay - maybe try some examples from their website. An error there would explain why no output is generated.
You can ignore any errors below line 5 of the above program output for now, as they are follow-on errors. In particular the question title isn't the 'real' error here, the real error is that no PDB files are generated on line 160.
By the way, I wouldn't use 10 500 as parameters for iteration number and number of models. It will take ages to run. I'd just use the defaults (omit the arguments), which we find gives good results.
Finally, I'll point out that the DMPfold repository has an issue tracker where you can raise any future issues. Though we may as well finish debugging this one here since we've started!

Ok, I cant ask more questions on this thread as stated by the person below. So, I will be asking next questions in another thread with title "DMPfold: is it suitable to get 5000 models to compare with rosetta".

thanks! I figured out the solution, it was the cns. I downloaded the binaries instead of compiling from source and it works!

Thanks for posting your solution! But please ask a new question for the follow-up question. This site is not a forum, so we don't do discussions, only questions and answers. You might want to take the [tour] or visit our [help] to understand more. In the mean time, I have removed the question from this answer, but please feel free to ask it as a new question.


text:
A program "PSICOV" "stopped" from a few days on Ubuntu 18.04?
I am running dmpfold, "seq2maps.sh" and from a few days the program is at the "PSICOV" stage:
psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

"ps aux" gives me this:
daniel     26693  490  2.5 2546588 2514320 ?     Rl   سبت20 25116:58 psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

Any idea on what could be the matter? Yes, I am runinng "rosetta" too on my 32 CPUs, so all my CPUs are covered to max 100%; howver, RAM is 93Gb and only 32 GB has been used. Thanks

In summary just be patient. psicov doesn't  parrallelise so a tenfold increase in CPU will not hurry it along.
Also psicov is not RAM intensive given the machine you are describing, so you haven't hit a RAM bottleneck.
I don't know the ps output format you are using I'm use to seeing %CPU and %MEM from ps -aux on Susie and OSX. You should also try,
top -c -u username

or for OSX top -u
You should see the psicov happily chuntering away consuming ~100% of a single process. There is a UNIX/Linux guru specialist on BioSE who could almost certainly provide better diagnostics.


text:
Can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?
I am trying to perform de novo structure prediction of a 1135 long protein sequence; can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?

For DMPfold you might like to look at Figure 5B of the paper, which implies that performance starts to drop off after ~600 residues. Only proteins up to 500 residues were used in the training set, for example.
A protein of that length is almost certainly multi-domain. I would consider doing domain prediction; we have a tool DomPred, but many others are available. Then run the sequences of the domains individually and recombine them later.
I don't know about I-Tasser's ability to do this, perhaps see if they mention it in their papers.

Generally, but not always the I-Tasser algorithm is not in the same league as a covariant mutation-based analysis (psicov) and don't think many people would disagree with this. I don't do I-Tasser but I'm aware it is the algorithm of choice for rough and ready predictions.
Psicov based models are essentially de novo predictions, based on amino acid mutation patterns, which I assume are refined against known homologous structure in DMFold (inaccuracies can creep through). I suspect the starting point of I-Tasser is a known homologue and it's associated known structure.
Having said that you need alot of sequence data to make psicov covariance work (>8000 aligned sequences) and (strictly in my opinion) it will also depend on the speed of mutation. Higher eukaryote  gene alignments could require much higher numbers of taxa in the alignment. The covariance value will give you the strength of the signal and values >0.7 should be used. I personally consider>0.5 to 0.7 borderline, which is fair enough because those are the generic thresholds in statistics. I don't know if you can access the underlying covariance values in DMFold, or moreover set the threshold.
So I-Tasser could outperform DMFold in given situations, but the underlying output of covariance based approaches should make it clear how well the approach has performed. Obviously if it performs well I'd go with that.
I would follow@jgreener 's advice on the follow up program.


text:
MiXCR: only create a single export file for all clonotypes
I am using following command from MiXCR to both align, assemble and export my input files:
mixcr analyze amplicon \
    -s hsa \
    --starting-material rna \
    --5-end no-v-primers \
    --3-end c-primers \
    --adapters adapters-present \
    --export "-aaFeature CDR1 -nFeature CDR1 -aaFeature CDR2 \
              -nFeature CDR2 -aaFeature CDR3 -nFeature CDR3 \
              -aaFeature VGene -cGenes -count" \
    --assemble "-OassemblingFeatures=[CDR1,CDR2,CDR3]" \
    INPUT_FILE.fasta analysis_name

Unfortunately, this creates multiple output file (analysis_name.clonotypes.ALL.txt, analysis_name.clonotypes.IGH.txt, etc.), which stands in contrast to the output of the mixcr exportAlignments command.
How can I export a single txt file instead, while still using the mixcr analyze command?

This is from the MiXCR documentation:
Under the hood analyze amplicon is equivalent to the execution of the following
MiXCR actions options:
> mixcr align -s hs -p kAligner2 input_R1.fastq input_R2.fastq alignments.vdjca
... Building alignments
> mixcr assemble alignments.vdjca clones.clns
... Assembling clones
> mixcr exportClones --chains IGH clones.clns clones.txt
... Exporting clones to tab-delimited file

So when using analyze amplicon, you are using mixcr exportClones under the hood and not exportAlignments, hence the difference in output.
I believe what you are experiencing results from a bug, please see this issue. If this is the case, all these files should be the same according to the bug report. Fortunately the developers have already addressed the bug in version 3.0.3.

As I was already running version 3.0.9 of mixcr, the bug you found (which corresponds to my observation) should have been fixed. However, although running mixcr exportClones prefix.clna clones.txt returns only one file as expected, the same is not true for me when using the analyze amplicon command. So, even though those commands are supposed to be equivalent, their output suggests otherwise for me.

Currently in MiXCR SERIES 4 you can use --dont-split-files  parameter here.
Please check our new documentation portal:
https://docs.milaboratories.com/mixcr/reference/mixcr-export/#clonotype-tables
Also now we support presets that allow to perform analysis for a lot of standard protocols using a one line command.

After finding the correct entry in the MiXCR documentation (thanks to haci), I learned that the parameter --receptor-type is responsible for which sequences will be returned.
The default for this parameter when using mixcr analyze is --receptor-type xcr, and stands for 'all chains', thus returning such a vast output. Changing this parameter to e.g. BCR solved my issue.

Where can I find the default value for the --receptor-type parameter? I cannot find it in the manual, can you please point out the relevant section?

I found it in the analyze section of the MiXCR docs, under the optional parameters.


text:
Why does this modification site file contain duplicate sites and multiple related proteins?
The data below is coming from a MaxQuant post-translational protein modification (PTM) site file.
Why does it contain duplicate sites (yellow) and proteins (cyan)?

Source data
http://ftp.ebi.ac.uk/pride-archive/2020/12/PXD020517/MaxQuant_BA39_SOLUBLE.zip

There are multiple protein isoforms that the gene may encode. The hyphen-number after the Uniprot ID represents the isoform —in a Uniprot you find the detail of these at the bottom.
The canonical isoform in Uniprot is the most likely one to be biologically relevant, while in NCBI and Ensembl it is the longest.
In a limited amount of case multiple isoforms are biologically relevant.
For the sake of sanity most variant annotation tools have a canonical argument to use the canonical.

According to a proteomics researcher at the NIH:

proteins column contains uniprot IDs.
position within proteins column are the starting position of each ID.
position is the leading protein


text:
Is loss/gain of function reflected in RNA-seq transcript counts?
Do LoF/GoF transcripts count toward the RNA-seq TPM count? Or would these LoF/GoF transcripts only be detected by isoform quantification?

Hi Kermit. Can you be more specific about what you mean by "loss-of-function/gain-of-function transcripts"? What scenario are you considering? Are there mutations within the coding regions of the gene, or are there transcript isoforms being expressed that change the function of that RNA? It's tough to know what you're asking about just from what you've written. Feel free to edit the question to add more clarifying information.

@JamesHawley The DNA mutation causes a big enough change that the final protein either no longer works or it takes on a different function

Short answer
Since the effects of any given mutation is so context-dependent (and deciding which effects matter will change between researchers), the process of "variant annotation" is typically separate from RNA-seq mapping/quantification.
If you use a common RNA-seq mapping/quantification tool like STAR, HISAT2, or Kallisto, these counts will combine transcripts with and without the mutation of interest.
This is because for the purpose of mapping, mutations are no different than sequencing errors, and when you're mapping you don't actually know what is what.
Long answer
You might not be able to quantify GoF/LoF at the isoform level, either.
Again, this is because the effects of mutations are so context-dependent.
If you want to know how many transcripts do or don't contain a specific mutation, you'll have to work backwards from that mutation and look at the RNA-seq reads themselves.
You can use tools like Rsamtools, htseq-count, or even IGV for this.
If you want to know whether a given protein's function is affected, you have to consider all the types of mutations than can lead to that protein being affected (in your biological scenario of interest, of course).

Some proteins may lose their binding affinity if a single amino acid is changed, in which case you can look for non-synonymous mutations in the RNA-seq reads.
Some genes may undergo alternative splicing and produce different proteins because of some mutations.
You could find this with a differential isoform expression analysis.
Some genes could be deleted on one allele, leading to a lower expression of that protein.
Biologically, that protein is still expressed, but having less of it makes it less functional in general.
You might be able to find that with copy number estimation or allele-specific expression analysis.

If you're looking at non-coding RNA, what mutations affect their function will be different from what you'd look for if you care about protein function.
All in all, any investigation about what a GoF/LoF mutation looks like and how to go about finding it will depend on you biological scenario.
And you may not find it by just looking at the TPM counts of a gene or its isoforms.


text:
blast nt vs Genbank
What is the difference between these two databases?
ftp://ftp.ncbi.nlm.nih.gov/blast/db/
 and
ftp://ftp.ncbi.nih.gov/genbank/
I understand that BLAST database has both nt for nucleotide and np for protein. But genbank is also for genomic database. So what is the difference between nt and Genbank databases?
Thanks

They are the same data, but the interegation of the database is very different. Blast uses a similarity algorithm but Genbank uses record referencing. It is better to limit the quesiton towards a specific biological problem for a more detailed answer

I think the BLAST folder contains pre-formatted databases for BLAST whereas the GenBank folder contains sequences.


text:
Which one is a more convenient assembly?
I have developed a software for de novo genome assembly. Its performance varies gradually according to how much data you employ. At initial stages it often produces contigs that look like that when aligned to a reference genome:

Initial-stage s-aligner assembly of SRR8357422
Genome fraction:98,8% -
Largest Alignment: 61.544 -
NGA50: 59.409 -
LGA50: 2
If I add more data, the contigs end up being larger, but that obviously takes more time. My question is, how does this result compare to, for example, one like that?

SPAdes assembly of SRR8357422
Genome fraction:93,46% -
Largest Alignment: 128.391 -
NGA50: 128.391 -
LGA50: 1
Note how you can construct a contig covering almost the full reference genome assembling green contigs in the first case. Therefore, a blast search of genes will likely let you find more genes given that it covers more genome fraction.
Apart from determining what option is more convenient, I would also like to know how much more convenient is to have both results. Or in which cases it's more convenient.
Last (but not least!), is there a way to generate automatically larger contigs from these overlapping contigs in assembly 1? I have tried using miniasm with no success. Even minimap2 seems not to find overlaps well enough. I am developing my own algorithm also for that but would like to know if there is another way.
[EDITED on 17/12 to add another example with real data for both assemblies]

Hi, could you please clarify your question a bit more; it's unclear to me. Asking "which one is a more convenient assembly?" suggests to me that there should be at least two comparisons, but you've only provided a figure for one. There is also additional information needed about the assembly and reference genome (i.e. are you doing kmer-based assembly; is the genome highly-repetitive; does it have Ns in the sequence).

Hi @gringer, it's not a question for a single case. It's something that happens in different assemblies in different circumstances. It's a generic question. The example in the image is real and corresponds to a run in PRJEB32127 but I have not assembled it with another software yet. I mean, in case I obtained an assemblie with a single contig but lower coverage (for example a contig 100k long) , which one of the two assemblies would be more useful.

Maybe it's more clear with a more extreme case. Suppose it's a bacterial assembly (~2Mbp) and SPAdes get a contig with 3Mbp bit just covering 90% of the reference (the rest is chimeric) while the otrher assembler gets lots of contigs around 200k bp long but covering 99%. I'd like to know how important is to get larger contigs or if coverage and percentage of genes covered is more important.

I will lsee if I find a real case with real data...

I guess it's "genome fraction" where I said "coverage"...

If it's a general question about assembly algorithms, it won't really fit well into the StackExchange Question/Answer format. Bioinformatics Stack Exchange works best for single, specific issues. I recommend posting your question on the bioinformatics subreddit.

Yes, makes sense. I will post it there. I will see if I find an example clear enough too and edit the answer or write a new one.

I asked that on Reddit. I didn't get an answer to all my questions but I got an answer for how to assemble the overlapping contigs.
If I use Flye I can assemble these and get larger contigs, better NG50, and higher genome fraction.

s-aligner + Flye assembly of SRR8357422 Genome fraction:95,95% - Largest Alignment: 134.103 - NGA50: 134.103 - LGA50: 1
Given that the data in the first assembly apart from being more useful itself in my opinion, is objectively more suitable to generate a better result in all senses with a simple additional step that takes a few seconds, my only possible conclusion is that it is overall a better result and therefore this sort of result will often be preferable.
Of course, there are other aspects to consider, like error rate, that are not considered here as it wasn't part of the question (mainly because I think Quast is not adequate for calculating it). But BTW, mismatches per 100kbp are also lower.
[EDIT] I want to add that this is not granted that in all cases Flye will generate a clearly better assembly.


text:
How to perform DE analysis for each sample
I am new to R and biocondunctor. I have the normalized expression values for 20 samples for a disease and for 10 controls. I wanted to get the differential expressed values for each sample with all 10 controls.
How could I do that with limma.
Any help will be appreciated.
Edit: What I want to achieve is single-sample t-test between the control population and the single sample values to get DE values

You need raw count data, and let limma normalize it (with e.g., voom). Read the manual, it is pretty well written.

@benn But I have the normalized data already, I want to get single expression values

Try to get raw counts.

@benn Ok , and then, I should normalize with voom, after which how should I get DE values for each sample separately?

Please go over the limma manual (https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf)  and edit your question with the potential problems you will encounter. Be careful to use the right section given the technology used (microarray or RNA-seq).

@haci Hey, I have added the microarray tag, I am sorry for the confusion

Okay I assumed it was RNAseq... Just read the manual...

@benn I am sorry if I have confused you, I tried to find this in the manual but I could not find any instance of comparing single sample vs a population

If you do statistics you compare groups, not separate samples. Separate samples do not have information about variance which is needed for statistics.

@benn But what I am doing is One Sample t Test. It may not have meaning as such, but I wanted to try it out with DE Values

Treat your disease samples as individual groups and then follow the normal routine in limma to use contrasts to compare two group (the control group versus individual disease samples). Note that you need to think long and hard about what these results then mean. Normally we compare groups because then the results should generalize to other samples. That will not be the case when you look at individual disease samples, so I suggest you do some clustering of the samples to see if you perhaps have any coherent clustering of the disease samples into subsets that could then be used as more useful groups in DE testing.

Thanks for the help!


text:
co-localisation error code using co-loc
Can anyone offer advice on what the below error code relates to. I am using co-loc R package with 2 GWAS datasets. Using coloc.abf function with pvalues, MAF, sample size and case/control cohort. Thank you
  results <- coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_data_subset$sample_size, type="cc",s=0.327,MAF=1_data_subset[,12]),   dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=0.4319,MAF=1_subset[,12], type="cc"),p12=1e-06)

error code: Error in data.frame(V, z, r, lABF) :
  arguments imply differing number of rows: 1994, 1992
In addition: Warning messages:
1: In 2 * N * f :
  longer object length is not a multiple of shorter object length
2: In 2 * N * f * (1 - f) :
  longer object length is not a multiple of shorter object length
3: In r * z^2 :
  longer object length is not a multiple of shorter object length

Hi, I can't replicate your error message. Can you double check that the line of code you inputted is the same as that which gave the error message.

The only thing different is that I put the code in a loop to make code quicker if(length(tmp)>0){
results = coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_subset$sample_size, type="cc",s=0.327,MAF=1_subset[,12]),
dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=CCRATIO,MAF=1_subset[,12], type="cc"),p12=1e-06)
}

In your dataset2 = list(), the MAF argument should be MAF=2_subset[,12] instead of MAF=1_subset[,12]. Since your two datasets have different numbers of lines, this is what causes the error.
I am not sure whether this is a typo in your code, or it is because you do not have MAF information in the second dataset. What I like to do before COLOC is to filter two datasets so they have the same number of lines, and then feed those filtered datasets to the coloc function.
(Also just out of curiosity, is 1_subset a data frame for your GWAS results? I think it is not recommended to have variable names starting with a number.)

Thank you for your reply. That is really kind. I dont have MAF for dataset two so used dataset 1 for both. I will try the filter function and then maybe input MAF detail first onto the second dataset using 1st dataset. I did change the names, didnt want to use the disease names on forum. Thank you for double checking! Have a safe day, kind regards!


text:
Parsing specific features from Genbank by label?
I'm trying to parse a Genbank file to find a specific feature. I can pull it out if I know the feature type (e.g. repeat_region). E.g. if I'm looking for this feature:
 repeat_region   5623..5756
                 /label=5' ITR
                 /note="5' ITR"

I know that I can find it using:
for feature in reference.features:
    if feature.type == "repeat_region":
        print(feature.location)

But I don't trust that it will always be a repeat_region. Instead, I'd like to look for it by label (5' ITR). I can seem to find a way to parse that from the feature object. Any suggestions?

The /label=5'ITR is called a qualifier. You can look for those qualifiers:
for feature in reference.features:
    for k, v in feature.qualifiers.items():
        if k == "label" and v == ["5' ITR"]:  # or "5' ITR" in v
            print(feature.location)


text:
match two files and check if position of chrom is within exon region
I have two files, file1 and file2. I want to match file1 chrom with file2 chrom, and if there is a match then I want to check file2 if the position of the match chrom1 is within the exon region that is in file1. i want to do this with the whole file and print out rows that match and the position is within the exonic regions.
For instance, file1 has a chrom 1 and POS 736689 which match chrom 1 in file2, now I take POS 736689 from file1 and check if it is within the exon region 29553,30563,30975,   30039,30667,31097,... in file2.
If it is within the region it will print out the file1 row that matches the pattern.
The exon regions can be read as:
29553 to 30039
30563 to 30667
30975 to 31097

file1 look like this:
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO
1   736689  T   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=6
1   758663  CGCC    TTAG    0/0 0/1 0|0 0|1 0|1 0|0 0|1 AC=2;AN=2184
1   762005  G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=13
1   771950  C   T   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=110;AN=2184
1   772172  C   T   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=8
1   790590  C   T   0/1 0/1 0|0 0/1 0/1 0|0 1|1 AN=2184;AC=84

file2 look like this
chrom   Start   End
chr1    29553,30563,30975,  30039,30667,31097,
chr1    53048,54829,    53067,54936,
chr1    69090,  70008,
chr1    160445,161313,  160690,161525,
chr1    324755,327551,  326514,328453,
chr1    367639, 368634,
chr1    459655,461153,461750,   459992,461288,461954,
chr1    523008,529614,  523081,530148,

Please edit your question so that your data is properly formatted (like you did withh file headers) and provide some easily parsable example data. I have posted an answer but not all people would (understandably) spend time cleaning data before suggesting a solution.

Despite an awk solution is favored by the OP, here is an R solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward awk solution, at least for me).
res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84

File 1 looks like a vcf file? Maybe you can read it in using VariantAnnotation package from Bioconductor, and do overlap with the second file. To make it easier, one can convert the second file to a GRanges object in R.

OP seems to add awk tags to all their posts even though multiple people have mentioned that R would be a better way to merge/lookup files.

It's a good answer by @haci nevertheless

Make a sorted BED file called exons.bed and run BEDOPS bedmap to map variants to exons:
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed <(vcf2bed < variants.vcf) > answer.bed

To make a sorted BED file from your exons, I'd suggest using sort-bed.
I would also suggest using awk to make sure that the chromosome naming scheme is consistent between exons and variants, e.g.:
$ vcf2bed < variants.vcf | awk '{ print 'chr$0' }' > variants.bed

Then use that result, instead, for mapping:
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed variants.bed > answer.bed

If chromosome names are not consistent, mapping or other set operations with BEDOPS or other toolkits will not work.

BEDTools already has a function to take a vcf and bedfile as input and look for intersections.  Most people would get the exon coordinates into bed format, and use that.


text:
How to check what software has been used to generate bam file?
I have several bam files from unknown origin. I'd like to know what software was used to obtain those bam files, e.g. bwa-mem. Is it possible to somehow check it?


text:
Finding bioperl on cpanmetadb failed
I failed to install bioperl on debian in the following way:
$ sudo cpanm bioperl
! Finding bioperl on cpanmetadb failed.
! Finding bioperl () on mirror http://www.cpan.org failed.
! Couldn't find module or a distribution bioperl

Thank you in advance,

conda install -c bioconda perl-bioperl  instead of using CPAN


text:
Mapping statistics from bam file using bbtools and sambamba
Below are the statistics for RNA-seq mapped and unmapped paired-end reads to rice genome using reformat.sh from bbtools on bam files. It gives 77% mapped and 5% unmapped, what about the remaining 18% reads? How I can get information for all reads?
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.mapped.bam, mappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.mapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 52020538 reads (77.10%) 7208830835 bases (77.41%)

Time:                         171.212 seconds.
Reads Processed:      67471k 394.08k reads/sec
Bases Processed:       9312m 54.39m bases/sec

For unmapped read;
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.unmapped.bam, unmappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.unmapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 3435326 reads (5.09%) 465456267 bases (5.00%)

Time:                         142.978 seconds.
Reads Processed:      67471k 471.90k reads/sec
Bases Processed:       9312m 65.13m bases/sec

With sambamba; Can I say 94.91% reads are mapped and 5.09% are unmapped?
67471075 + 0 in total (QC-passed reads + QC-failed reads)
12015211 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
64035749 + 0 mapped (94.91%:N/A)
55455864 + 0 paired in sequencing
27727932 + 0 read1
27727932 + 0 read2
51006046 + 0 properly paired (91.98%:N/A)
51419222 + 0 with itself and mate mapped
601316 + 0 singletons (1.08%:N/A)
263816 + 0 with mate mapped to a different chr
237111 + 0 with mate mapped to a different chr (mapQ>=5)

Try samtools flagstat or sambamba flagstat to get similar type of information about your bam file. It may be easier to interpret. Both of these tools can be installed using conda if needed.
EDIT:
Yes, the results of sambamba flagstat you showed in your edited post indicate exactly what you said: 94.91% mapped and 100% - 94.91% = 5.09% unmapped:
64035749 + 0 mapped (94.91%:N/A)

Hi Timur, Thanks, I have added the results from sambamba as well, I think it gives the accurate result. Could you please respond to query related to sambamba result?

@bioinfonext Edited the answer.

More of a comment than an answer... In my opinion the terminology is a bit confusing since the term reads actually refers to a more generic records.
Look at this test file:
samtools view test.sam | cut -f 1-4
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529094
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:403:C6F0HANXX:8:2215:4792:20065 4   chr7    5529871

I would say that here there are 5 records (of which 4 are alignments) but only 2 reads.
samtools flagstat reports 80% mapping rate but maybe I would say this is 50%.
samtools flagstat test.sam
5 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
4 + 0 mapped (80.00% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)


text:
insert size pre and post trimming
I have a problem here with my rna seq data:
Sequencing details:
rRNA was removed, followed by cDNA preparation and generation of stranded libraries using the TruSeq Stranded Total RNA Sample Prep Kit. Sequencing performed on the HiSeq2500 platform (Illumina) to generate 2 × 125 bp paired-end reads
Alignment and preprocessing
reads were aligned using tophat2(std aligner in pipeline at core (--library-type fr-firststrand), unfortunately the original unaligned files were purged and i only have access to this aligned file which i converted to fastq using following steps

samtools sort -n sample.bam -o sample_sorted.bam
bedtools bamtofastq -i sample_sorted.bam -fq sample_1.fq -fq2 sample_2.fq (get a lot of mate skipping errors here)
check for these reads for adapters and trim them and again use hisat2 (with --rna-strandness RF)

trimmed reads
Reads were trimmed by passing parameters --adapters adapters.file --adapter-trim-end RIGHT --length-dist --threads 12 --adapter-min-overlap 7 --max-uncalled 250 --min-read-length 25 -- to FLEXBAR version 2.4
adapters for all samples using multiqc

fastqc sequence length pre trimming

fastqc sequence length post trimming

rmats error

Incorrect readLength. sample.bam has a read length of 114, while readLength param is 125

picard error on aligned bam from hisat2)

ProcessExecutor   2 "Not creating insert size PDF as there are duplicated header names: All_Reads"

snapshot of picard insertsize metrics file

Questions

Is my approach flawed at some step?
why do i get different read lengths error for each sample when i process it through rmats when i specify len as 125?
why is no histogram being produced when using picard insertsize metrics?
is normal to see the graph shift to the right after trimming?

insertsize for all samples (data from picard collectinsertsize compiled using multiqc)

It's not an unreasonable approach. rMATs is rather picky about its input, but it seems you've noticed that.
rMATs can't handle trimmed reads. It also can't handle soft-clipped reads. You might as well not trim. You'll get a lower alignment rate, but your only other choice would be to either exclude the trimmed reads or trim everything down to the same length.
You have a mixture of RF and FR alignments and the plotting script can't handle that.
Assuming you mean the read length distribution, it's actually shifted to the left, not the right. You had a single length before and a long left tail after trimming.

How does one get a mixture of RF and FR, when the library was FR ? doesn't that mean the alignment is doubt full as well?

@novicebioinforesearcher That's a very good question and not one that I have a good answer for.

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


text:
How can I export a full alignment from IGV as an image?
The IGV browser lets you export an alignment as an image (File => Save Image). However, this image only contains those reads that fit in the viewing window:

As you can see in the image above by looking at the scroll bar on the right, this is only a small subset of aligned reads. Is there any way to export an image of the aligned reads that includes all reads and not only those currently shown in the window?

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue. To run a script, first make a script file like so with 1 command per line: snapshotDirectory ./ goto chr2:28582094-28582790 snapshot chr2_28582094_28582790.png goto chr4:8203487-8203815 snapshot chr4_8203487_8203815.png Run the script using Tools | Run batch script menu.

@TimurShtatland sorry, I have no idea what that is supposed to mean. That looks like a set of scripting commands and I have a vague recollection of having seen something like that for IGV in the past. Could you please post an answer explaining what you are suggesting?

Oops, pasted the wrong thing. Fixed this. IMO, this is not an answer, but might help partially.

@TimurShtatland yes, thanks. I'd posted my comment before your edit :). If the scripting approach does give more reads, then it's at least a partial solution and deserves an answer. Could you flesh your comment into an answer? If you'd rather not, I can post one as soon as I've tested it, but you deserve the rep.

@TimurShtatland well, this is embarrassing. Your comment made me find an old script I had which does exactly this! I had completely forgotten! Please post an answer, you should get the credit, and then I'll just edit it to add what I did to get everything (essentially, you can set maxPanelHeight to get more reads).

Thank you, done. Please feel free to edit. I don't deserve the full credit obviously, so I owe you half. :)

@TimurShtatland thanks, but that's fine. The credit belongs to Pierre whose post I'd managed to find last time I needed this and then promptly forgot about until your comment reminded me :). And I already have all the editing privileges rep can buy here, so you need it more!

Let us continue this discussion in chat.

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue by itself.
To run a script, first make a script file like so with 1 command per line:
snapshotDirectory ./
goto chr2:28582094-28582790
maxPanelHeight 10000
snapshot chr2_28582094_28582790.png

Run the script using Tools | Run batch script menu. The importat trick is the maxPanelHeight 10000 which is what makes this include more reads than the default snapshot.
You can also convert this into a little shell script:
#!/bin/bash
mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
goto chr2:28582094-28582790
sleep 2
echo "maxPanelHeight 10000"
echo "snapshot chr2_28582094_28582790.png"

Then, just open IGV, load the relevant BAM file and run the script like this to connect to IGV via telnet (credit to Pierre Lindenbaum):
sh foo.sh |  telnet 127.0.0.1 60151


text:
How do you generate read-length vs read-quality plot for long-read sequencing data (e.g., MinION)?
How do you generate read-length vs read-quality plot (heat map with histograms in the margin) for long-read sequencing data from the Oxford Nanopore Technologies (ONT) MinION? The MinKNOW software from ONT provides a plot like this during base calling.
This would also be very helpful for PacBio data.

I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:

A fastq file (optionally compressed)
A bam file
The sequencing_summary.txt file generated by albacore

Using optional flags you can:

Log transform the read lengths
Use aligned reads rather than sequenced reads
Downsample the reads
Set a maximum read length

I've added an example below, plotting log transformed read length vs average read quality (using a kernel density estimate). More examples can be found in the gallery on my blog.
I welcome all feedback and suggestions!

It's important to always consider read length and quality jointly with high-error read data, and current long-read technologies (e.g., MinION and PacBio) have high error rates. Considering read length and quality jointly will help you determine how successful the run was, how many reads were 'high quality', whether the longer reads are 'real' (or just pore noise), etc.
I've had a recent spike of interest in similar plots and came across a project called pauvre (french for 'poor', play on 'pore') through the Oxford Nanopore Technologies (ONT) community that I think is even better than MinKNOW's base calling plot.  Plus, you can generate these plots from a fastq file when ever you want to, unlike with MinKNOW.
[Note: I am not the original author, but I am now contributing because I liked (and needed) it.]

Pauvre will also report useful statistics:
fastq stats for fastq_runid_bb8b8ddedb22bdd6802b2bfa2b4e424c92c30d28_0.fastq
numReads: 2164829
numBasepairs: 4970615217
meanLen: 2296.077527139557
medianLen: 1495.0
minLen: 5
maxLen: 392031
N50: 3450
L50: 402786

                               Basepairs >= bin by mean PHRED and length
minLen          Q0          Q5         Q10         Q15       Q17.5         Q20      Q21.5  Q25  Q25.5  Q30
     0  4970615217  4970611559  4835461771  3889995868  2900103275  1087779109  165637656  429      0    0
 50000    11531044    11531044      270324      160128       50729       50729          0    0      0    0
100000     6260554     6260554           0           0           0           0          0    0      0    0
150000     3504240     3504240           0           0           0           0          0    0      0    0
200000     2501101     2501101           0           0           0           0          0    0      0    0
250000     1609592     1609592           0           0           0           0          0    0      0    0
300000     1033423     1033423           0           0           0           0          0    0      0    0
350000      392031      392031           0           0           0           0          0    0      0    0

                    Number of reads >= bin by mean Phred+Len
minLen       Q0       Q5      Q10      Q15    Q17.5     Q20  Q21.5  Q25  Q25.5  Q30
     0  2164829  2164605  2083436  1626706  1183812  435687  77341    1      0    0
 50000      109      109        5        3        1       1      0    0      0    0
100000       36       36        0        0        0       0      0    0      0    0
150000       15       15        0        0        0       0      0    0      0    0
200000        9        9        0        0        0       0      0    0      0    0
250000        5        5        0        0        0       0      0    0      0    0
300000        3        3        0        0        0       0      0    0      0    0
350000        1        1        0        0        0       0      0    0      0    0

These plots and stats would be equally useful with PacBio, but that's not super easy (though it is possible) with current raw output from the Sequel sequencer: Which quality score encoding does PacBio use?
Pauvre currently uses Biopython to parse the fastq and matplotlib for the actual plot, and will let you choose the output image format (e.g., .png, .pdf, etc.). You can also choose whether the background is transparent or white (for .png output).
The parser is currently super slow because it's using SeqIO.parse, but we're changing parsers to speed that up. We're also adding some extra features (e.g., choose whether to include y-axes in margin histograms, print some stats directly to the plot for documentation, etc.)
Purple is currently the only color choice (which I personally love), but adding options to change that will be super easy.

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


text:
Get raw genome sequence
I'm currently working on testing some matching algorithms for strings. I would like to do some tests on raw genomic data as I expect different results from random strings given the lower entropy. How can I get a whole human genome sequence? (just the raw basis)
I have tried websites like https://www.ncbi.nlm.nih.gov/ it seems very confusing and full of unnecessary information given that I just want to have the raw basis without annotations.

Its potentially an interesting project. I'm aware of comparable projects, but ... your question and investigation is ill-defined. Simply if you download 'fasta' format for any given sequence you will get the description on one line beginning with '>' and the sequence on the line below. Genbank is a good place to go

@Michael I'd be happy to have references to the comparable projects, thanks!

Just a small comment, "raw" usually in genomics means unprocessed sequencing reads: https://en.wikipedia.org/wiki/Read_(biology), the sequence you seem to be after is "genome assembly".

You want FASTA files from here.
If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).
See the README for more info.


text:
How can I extract information from .sam files?
I have 10 .sam files after my bowtie2 alignment on ten single-pair sequences. I would like to build a graph based on that output data, however I don't see any other solution but just writing the alignment percentage manually after each iteration.
Is there a tool that can help me to build a .csv file so I can manipulate the output data more easily?

You should never have SAM files: they’re not useful, not even as intermediate files! You should always directly generate BAM files. They’re smaller, more efficient, and allow more features (e.g. indexing for even faster access and statistics). For example, you would usually call bowtie2 as follows: bowtie2 ‹args› | samtools sort --output-fmt-option level=0 | samtools view -b -o sorted.bam (this immediately sorts the BAM file, which requires a lot of RAM; if this isn’t suitable for you, omit the intermediate step; but you normally want to have sorted BAMs).

Use samtools flagstat with option -O tsv:

Using -O tsv selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names.

SEE ALSO:
What Does Samtools Flagstat Results Mean?
INSTALLATION:
Some old versions of samtools do not have option -O tsv. To install a recent versions of samtools, use either:
conda (macOS, Linux) (but see also: Tutorial: installing samtools with conda - overcoming two common errors):
conda create --name samtools samtools

or brew (macOS):
brew install samtools

Thank you, it worked however just in a "default" combination, e.g. samtools flagstat input.sam           If I try  samtools flagstat -O tsv input.sam, it throws me an error  invalid option -- 'O'

@DmitriiTrubetskoy See updated answer.

thanks, however my  samtools is already the newest version (1.7-1). (info from terminal). Just to be sure, the right syntax will be: samtools flagstat -O tsv input.sam ?

@DmitriiTrubetskoy Exactly, and note that 1.7* is older than 1.11, which is the most recent version available through both conda and brew (see links in my answer). I also confirmed that this option is present in samtools 1.10 that I installed via brew. So try updating/installing the new samtools version if you really need this option.

Thanks a lot. I just wonder how can I roll back to 1.10? Is there a way to do it?

@DmitriiTrubetskoy Don't roll back, because 1.7 is fairly old. I am suggesting updating from 1.7 to the most recent version, currently 1.11. Use conda or brew. To install conda (for mac or linux), follow directions here for miniconda: https://docs.conda.io/projects/conda/en/latest/user-guide/install/ . To install brew (for mac), follow the directions here: https://brew.sh/ . Then use the commands in my answer to install samtools (via conda or brew). This will install by default the most recent version. You are usually best off with the most recent version. Good luck!

Thank you. Now it worked, something happened...however I suppose I must have created some new file that would contain some output information - what kind of file that should be? I cannot see find anything in the same folder where input.sam file was.

@DmitriiTrubetskoy It writes into STDOUT a tsv file, 16 rows x 3 cols.

Great, now everything's ok. Much appreciated the tips and for being patient.


text:
pysam or piping samtools view to a python script
Is it faster to use the PySam package to run a python script on a bam
for read in samfile.fetch('chr1', 100, 120):
 print read

compared to using a pipe and reading line by line from stdin in a python script.
samtools view example.bam | python script.py.

for line in tqdm(sys.stdin, total=number_of_lines)

It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will not be the bottleneck in your processing, rather the python script itself will be.


text:
Good / recommended way to archive fastq and bam files?
We have a lot of Illumina sequenced exome data. Currently we are using spring for its great lossless compression, but we are looking if there is anything better (and most preferably opensource) which can let us compress our fastq files.
We also want to compress or reduce file size of BAM, but using CRAM in lossless mode didn't yield good results. The space savings were decent, roughly ~10% (IMO was expecting more), but the biggest hit came when using it with IGV. We have set up a http-server and stream bams to IGV (Helped us immensely as people wouldn't have to download and visualize bam files). With Cram, the loading times were significantly high, requiring users to commit more memory for IGV. Sometimes we would have to wait around a minute for reads to show up. With Bam it is near instantaneous.
Can't use lossy compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).
What are some good ways to compress fastq and bam files so we can maximize our storage?

EDIT: I am rewriting the answer in response to updates to the original question.
TL;DR: use CRAM
Background 1: quality binning and FASTQ compression
In the old days, base callers outputted base quality at full resolution – you could see quality from Q2 to Q40 in full range. As a result, quality strings were like semi-random strings and very difficult to compress. Later people gradually realized that keeping base quality in low resolution wouldn't affect downstream analysis. The Illumina basecaller started to output quality in 8 distinct values and later changed that to 4 bins. This change greatly simplified quality string and made them compressed better. For example, in old days, a 30X BAM would take ~100 GB. With quality binning, it would only take ~60 GB.
Background 2: GATK Base quality recalibration
In early 2010s, Illumina base quality was not calibrated well. GATK people introduced BQSR to correct that and observed noticeable improvement in SNP accuracy. Nonetheless, with improved Illumina base caller, their base quality became more accurate. Meanwhile, the world moved to 30X deep sequencing. The depth overwhelms slight inaccuracy in quality. I would say around 2015, BQSR was already unnecessary for data produced at the time.
Does it hurt to apply BQSR? Yes. First, BQSR introduces subtle biases towards the reference and towards known SNPs. Second, BQSR distorts the data. At least for some datasets, I observed that SNP accuracy dropped with variant quality after BQSR; I didn't observe this with raw quality. Third, BQSR is slow. Fourth, for new sequencers producing data at higher quality, BQSR is likely to decrease data quality. Last, related to the question, BQSR added another semi-random quality string and made compression even harder. Nowadays, running BQSR is a waste of resource for worse results. The official GATK best practice no longer uses BQSR according to their WDL file.
CRAM
These days a 30X human CRAM only takes ~15 GB (see this file). This is a huge contrast to ~100 GB BAM in early 2010s. OP only saw ~10% saving probably due to a) BQSR and/or b) old data with full quality resolution.
On encoding/decoding speed, CRAM was much slower than BAM. Not any more. The latest htslib implementation of CRAM is faster than BAM on encoding and only slightly slower on decoding. The poor performance of IGV on CRAM could be that the Java CRAM decoder is not as optimized.
It is true that CRAM is not as widely supported as BAM. However, all the other alternatives are much worse. Petagene said they had IGV-PG (PDF) for their format. That is not official IGV and I couldn't find more recent update beyond the 2019 press release. I don't see other viable options.
Note that the common practice is to keep all raw reads, mapped or not, in BAM or CRAM such that you can get raw reads back later. BAM/CRAM additionally keeps metadata like read group, sample name, run information etc and is actually more popular than FASTQ in large sequencing centers. Also note that you don't need to sort CRAM by coordinate. Unsorted CRAM is only a little larger than sorted CRAM.
CRAM and its competitors
The core CRAM developer, James Bonfield, is one of the most knowledgeable researchers on compression (and one of the best C programmers) in this field. He has done a lot of compression evaluation over the years. The conclusion is that on a fair benchmark, CRAM is comparable to the best tools so far in terms of compression ratio.
Petagene could compress better in the plot @terdon showed mostly because it has a special treatment of the OQ tag generated by GATK BQSR. It is a typical trick marketing people use to make their methods look better. With BQSR phased out, this plot is no longer relevant.
On commercial software
In general, I welcome commercial tools and think they are invaluable to users. I also have huge respect to DRAGEN developers. However, on FASTQ storage, I would strongly recommend against closed-source compressors. If those tools go under, you may lose your data. Not worth it.

Thanks for your input, but mainly for the link. WDL seemed too complicated for me. We also noticed that BQSR was not recommended by nf-core people as well, but none provided a concrete reason beyond "its not useful". As I've also commented on @gringer post, CRAM didn't work well when using it with IGV, as well as some downstream tools don't work with CRAM.

@KarthikNair It is interesting that you are criticizing CRAM for slowness in IGV, while you are compressing fastq with spring and you can't even load other alignment formats into IGV. What are you expecting? Also, you don't need to understand WDL. Just search/grep "BaseRecalibrator" or "BQSR" to see if it gets called.

@KarthikNair sadly, the GATK best practices, which were updated very recently, still recommend BQSR so it is really surprising that their recommended pipeline doesn't use it. That said, in my place of work, we analyze around 2k human samples a month, and we have found that BQSR often causes us to miss variants. I don't have statistics for this yet so it's anecdotal, but there is strong evidence that BQSR not only isn't particularly helpful today, but can be actively harmful. And not only with small panels, this also happens with whole exome sequencing samples.

@terdon That page is referencing their older pipelines, which do use BQSR. The WDL file in my link is two years old. It is the latest I can find. IMO, BQSR has been hurting for 10 years. Well, better late than never.

you focus on BAM but what would you suggest for FASTQ if not DRAGEN ORA reference based approach? Arguably long term of storage of FASTQ is more important since BAMs can be regenerated

@Chris_Rands See new edits. I see the value in reference-free algorithms such as spring and dsrc. But if you already require a reference, it is better to keep CRAM than ORA.

For what it's worth, petagene's compatibility layer lets both IGV and JBrowse work out of the box with peta-compressed bam files. You don't even need to use their IGV implementation, vanilla IGV works fine as long as you are using petagene's trick with the symlinks I mentioned in my answer.

@Chris_Rands You can align raw reads with STAR or any local aligner and then create a CRAM. I guess that will give you a smaller file than BAM.

Hay all, thanks for all your comments! I will echo @terdon and encourage people to share their findings w.r.t. to BQSR not being very helpful. We have been using MarkDuplicatesSpark and BQSR pipeline. I'm planning on dropping the BQSR steps and using samtools markdup -s instead of MarkDuplicatesSpark (Mainly for optimization). And yes, perhaps the lackluster size reduction could have been due to BQSR in data. We also have old NGS data sequenced from 2016 onwards, so some samples might use unbinned quality scores.

I'd recommend reconsidering CRAM. It's an open standard that is most likely to still be accessible many years in the future. Regarding your point about lossy compression...

Can't use lossy [CRAM] compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

By default, CRAM compression is not lossy.
CRAM does, however, have an optional archive settings mode (samtools view --output-fmt cram --output-fmt-option archive ...) which is a lossy compression, doing things like removing read names, removing additional accessory fields, and additional compression of quality scores. In all cases, the base sequence of the reads is preserved.
https://www.htslib.org/doc/samtools.html

An alternative to the good option of PetaGene that @terdon mentioned and also not open source is the DRAGEN ORA compression. This is loseless reference based compression and can be ~5 x smaller over fastq.gz, but it's only for fastq files and only those generated by Illumina machines (and will be easiest within the DRAGEN ecosystem).
Another point to consider is where you store the data. Even with 'standard' compression you might get more benefits moving archival data to a system like AWS S3 Glacier, depending on your current infastructure setup of course

This is exactly what I do. I archive FASTQ.GZ (with added content-md5 checksum) to S3. I also archive processed outputs to a different location that is less fault tolerant than S3 as those files can be recreated.

We did check this out, but we do not sequence the NGS data ourselves. We ask external labs to perform the sequencing and they share the fastq.gz files with us. What is the approximate pricing for this? I can't clearly identify where those details are!

The only free and open source tool I know that can help is zstd. Their github repository's README describes it as:

Zstandard, or zstd as short version, is a fast lossless compression algorithm, targeting real-time compression scenarios at zlib-level and better compression ratios. It's backed by a very fast entropy stage, provided by Huff0 and FSE library.

We read this blog post on using it to compress fastq files and were intrigued so we ran some tests. We couldn't reproduce the same level of compression that the article claims but that might be due to different input files or the way an index was/wasn't created.
Even in our tests though, we did see that it provides better compression than gzip. Here are the specific values for one of the tests, run using an R1 file of a human WGS run:

Uncompressed fastq
bgzip-compressed
zstd-compressed

179G
58G
47G

This was done using -z 19, to set the compression level to 19 (we found this to be the best tradeoff between size and time to compress). We also saw that using a human genome fasta file to create a compression index, and then using that when compressing human fastq files, actually resulted in worse performance (larger compressed file) than not using the index. But it is entirely possible that we just did it wrong, we didn't spend a lot of time on it.
So you might want to investigate zstd as well, maybe you can get better performance if you play around with the indexing. On the plus side, it is free
and open source. On the other hand, it won't do much with bam or cram files.

Thanks for the input. I'm pretty sure that with Spring, the output size might be around 20-30 GB, if not slightly lesser. Granted, you would have to compress / decompress the data, but the final file size would be pretty small. We have compressed WGS samples with spring, it goes from total 50GB of gz files to around 15-20 GB in final spring file.

Zstd (and Zpaq and any other general compression tool) deliver relatively bad performance on SAM/BAM files. This is because they are either LZ-based or otherwise view the text as being "one long text", whereas SAM/BAM is anything but that.
FYI, this opinion is based on fairly thorough testing.

You might want to consider PetaGene. This is a paid tool, and not open source, but it is something we use at work and it is really good! It can do lossless compression of gzip-compressed fastq files, bam files, and even cram files, and achieves really good size reduction (image from the PetaGene webpage):

In my own testing, using files of different sizes, I get an average size reduction of 62.1% (2.6 x smaller) for fastq.gz files (compared to the gzipped fastq file, not compared to the original, uncompressed fastq file) and an average reduction of 82.1% (5.6 x smaller) for bam files. I didn't test cram files.
If you have the budget for a paid service, these guys do good work! The greatest advantage, apart from the very significant size reduction, is that they also provide a system that allows you to use the compressed files natively. So you can feed the petagene-compressed bam file directly to bwa for example. They create special symlinks (foo.fastq.gz -> foo.fastq.petacompressed) and you can then use these symlinks and petagene handles the decompression on the fly.

The Wikipedia page on CRAM has a nice table with a few alternatives to CRAM which I reproduce below. I haven't tried any except PetaGene though so cannot comment on their performance:

Name
Key advantage over CRAM
Repo

Deez
None. An early approach to SAM compression, introduced 2014.
https://github.com/sfu-compbio/deez

Genozip
Significantly better compression, but slower random access
https://github.com/divonlan/genozip

GenomSys
Utilizes ISO-standard MPEG-G, but compression inferior to CRAM 3.1
Source code not publicly available

PetaSuite
Virtual filesystem interface
Source code not publicly available

We considered this software, but sadly we don't have the budget for this. And yes, on-the-fly decompression is a great feature to have, makes storing all data in this format extremely viable. Right now, we get decent compression with spring, but it needs lots of de/compression time, not to mention the CPU resources it guzzles during the process lol

@terdon genozip is quality work and is very impressive. Its full source code is also on github. If genozip goes under some day (I hope not!), you may create a fork and continue to use their format.

@user172818 ~terdon - can you share any further info on how PetaGene achieves compression and the critique mentioned of "special treatment of the OQ tag generated by GATK BQSR"? I know it's closed source but trying to understand conceptually how their compression works

Not I, @Chris_Rands, sorry. I'm sure user172818 will know far more, but I have no idea and wasn't even aware there was any special treatment before reading their answer. I know at least one person who used to be an engineer working for them is active on the site. If there's something that can be shared, maybe they will answer if you post a question.


text:
Nextflow Error: failed to read header from "-"
I am trying to run my nextflow pipeline, and have gotten this error:
samtools sort: failed to read header from "-"
I'm not sure why this error is occurring.
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
import groovy.io.*

params.t = 8

process joinFastq {
input:
path(ds_directory)

output:
tuple path("joined_${ds_directory.baseName}.fastq"), val("$ds_directory.baseName")

script:
"""
gunzip -c ${ds_directory}/*/fastq_pass/*.fastq.gz >
 joined_${ds_directory.baseName}.fastq
"""

}
 process minimap2 {
cpus 8
memory "20 GB"
input:
tuple path(fastq), val(baseName)

output:
tuple path("*.bam"), val(baseName)

  script:
   """
  minimap2 -ax splice -uf -k14 -t 8 ${params.genome} ${fastq} > mapped_${baseName}.sam
    samtools view -h -b mapped_${baseName}.sam > ${baseName}.bam
   samtools sort -o ${baseName}.bam > ${baseName}_sorted.bam
   """
   }

I have tried applying the -h flag to samtools sort in the script
   samtools sort -h -o ${baseName}.bam > ${baseName}_sorted.bam

and the same error is produced. I'm not sure why the header flag didn't solve this issue.

This is not a Nextflow error. The problem is that without an input file, samtools sort tries to read from stdin. Using a recent samtools, you can however coordinate sort the SAM and write a sorted BAM using:
samtools sort -o "${baseName}.bam" "mapped_${baseName}.sam"

You may have been intending to pipe the output to samtools sort, which would avoid writing large SAM files and is usually preferable. A minimal example might look like:
params.reads = './*.fastq.gz'
params.reference = './gencode.v43.transcripts.fa'

process minimap2 {

    tag { sample }

    cpus 8
    memory 20.GB

    input:
    tuple val(sample), path(fastq)
    path ref_fasta

    output:
    tuple val(sample), path("${sample}.bam")

    script:
    def cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    minimap2 \\
        -ax splice \\
        -uf \\
        -k14 \\
        -t ${cpus} \\
        "${ref_fasta}" \\
        "${fastq}" |
    samtools sort \\
        -o "${sample}.bam"
    """
}

workflow {

    reads = Channel.fromFilePairs( params.reads, size: 1 )
    reference = file( params.reference )

    minimap2( reads, reference )
}

Also, be sure to avoid using absolute paths in your workflow scripts. For your workflow to be portable, processes should only access files in the working directory. The example above avoids this, by using a value channel, to stage the reference FASTA file into the working directory.
In case anyone is wondering, the ternary is to ensure that we do not oversubscribe the node(s) that these jobs land on. Otherwise, each job could use (up to) 900% (800% + 100% = 900%) CPU.


text:
How to load and split my two individual datasets when integrating datasets in Seurat?
I want to do comparative analyses of my two individual single cell RNA-Seq datasets in Seurat. And I load and list my data as follows:
Control.data <- Read10X("~/Documents/...  Control_filtered_feature_bc_matrix")

Control <- CreateSeuratObject(counts  = Control.data, project = "Control_summary", min.cells = 3, min.features = 200)

Mut.data <- Read10X("~/Documents/... /Mut_filtered_feature_bc_matrix")

Mut <- CreateSeuratObject(counts  = Mut.data, project = "Mut_summary", min.cells = 3, min.features = 200)

Object<- list(Control,Mut)

names(Object) <- colnames(Object)

Then the protocol https://satijalab.org/seurat/v3.1/immune_alignment.html said, I need to use SplitObject function, but what do I split by?
Object_list <- SplitObject(Object, split.by = "?")

Can you please edit your post so as to include the link for the protocol that you have mentioned? Moreover, you have formatted your code as a "quote", please also change this to "code".

Seurat is expecting individual datasets to be normalized separately prior to data integration. In that respect lists of objects corresponding to different datasets are handy to manipulate each object/dataset individually.
In your case you do not need to do any splitting, you have an object corresponding to each of your datasets.
If you would have started with a single count matrix of all datasets then you would need to split your data based on batches/experimental condition/...


text:
awk find pattern
I am trying to find this pattern on a large file but the output doesn't look good. The output doesn't seem to reflect the pattern.
Father = 1|0 or 1/0

Mother = 110 or 1/0

Daughter 1 = 1|0 or 1/0

Daughter 2 = 1|1 or 1/1

Daughter 3 = 0|0 or 0/0

Son 1 = 1|0 or 1/0

Son 2 = 1|1 or 1/1

The content of the file1 looks like this
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3

SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   3833321 T   C   1/0 1/0 1|1 1/0 1/0 1|1 0|0 AN=2184;AC=801  0.366758242
1   9009352 A   G   0/1 0/1 1|1 0/1 0/1 1|1 0|0 AC=90;AN=2184   0.041208791
1   9328077 G   A   1/1 0/0 1|0 1|0 1|0 1|0 1|0 AC=1;AN=2184    0.000457875
1   10684412    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=2    0.000915751
1   11782120    G   A   1/0 0/0 0|0 0|0 1|0 1|0 0|0 AC=3;AN=2184    0.001373626
1   12034879    A   T   0/1 1/1 0|1 0|1 0|1 0|1 1|1 AC=3;AN=2184    0.001373626
1   16064465    C   A   1/0 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=4    0.001831502
1   16111170    G   A   1/0 1/0 1/0 0|0 0|0 1/0 1/0 AC=1;AN=2184    0.000457875
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   16258118    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=1    0.000457875
1   17966719    A   G   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=0;AN=2184    0
1   19751079    CG  C   0/1 0/1 1|1 0|0 0|0 0/1 0/1 AN=2184;AC=1    0.000457875
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
1   21060356    G   A   1/0 1/0 1/0 1/0 1/0 1|1 0|0 AN=2184;AC=541  0.247710623

My awk script look like this
awk '($8~/1\/1/ || $8~/1\\|1/) && ($9~/0\/0/ ||$9~/0\\|0/) && ($7~/^1\\|0/||$7~/1\/0/) && ($10~/1\\|0/||$7~/1\/0/) && ($10~/0\/0/ ||  $10~/0\\|0/) && ($5~/1\/0/|| $5~/1\\|0/) && ($11~/0\/1/||$11~/0\\|1/)&& ($6~/1\\|0/||$6~/1\/0/) && ($12 ~/1\/1/ ||12~/1\\|1/){print}' file1 >file2

What's wrong with the output? "Doesn't look good" isn't very informative. Please [edit] your question and highlight what is going wrong.

The output doesn't reflect the pattern.

Well yes, I imagine so. But since this is a really complex pattern, it would be very helpful if you could show the output you get and explain why it is wrong. I wrote you the script you're using, and I still don't see what the issue is.

Example output of what you got and what you expected would be helpful, without these one can only guess.
First of all your input does not have a row satisfying the criteria you have given, so I added two dummy rows, one with / and one with | used as separator.
Here is the code I used:
awk '($5~/1\|0/ || $5~/1\/0/) && ($6~/1\|0/ || $6~/1\/0/) && ($7~/1\|0/ || $7~/1\/0/) && ($8~/1\|1/ || $8~/1\/1/) && ($9~/0\|0/ || $9~/0\/0/) && ($10~/1\|0/ || $10~/1\/0/) && ($11~/1\|1/ || $11~/1\/1/) {print}' test.txt

And its output (input is yours + 2 dummy rows):
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy

The code is basically yours, except for:

I escaped | with one \ instead of two
I changed this part: ($10~/1\\|0/||$7~/1\/0/). I think $7 should have been $10
I removed a ^ from one of your conditions, I believe they all should start with ^ though, per @terdon's answer on a related question

I believe the code above does what you want, but when I run code above with two backslashes to escape I get the following output and could not figure out why!
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy


text:
Convert PAF format to SAM/BAM format
I have a bunch of PAF files resulting from the alignments of fastq files on a reference genome with minimap2. I would like to convert them into SAM/BAM format so I can use samtools mpileup to recover the read coverage at filtered positions
So far I've looked at paftools.js and found sam2paf but no paf2sam
I am open to suggestions but I'd like to avoid mapping the reads again.
Edit
I have opened an issue on minimap2's github
Answer from Heng Li :

You need to rerun alignment. No tools exist to convert paf to sam.
why there is no paf2sam in paftools.js

It is technically possible, but someone needs to write the code.

I'm not sure the question is answerable in its current format, there doesn't appear to be enough background information and this often results in substancial misunderstanding between OP and the reply. Lets see what the other think.

Thanks for the feedback. So, as far I as know mpileup can take an alignment file and a list of genomic positions (eg BED) and return statistics for each positions, such as the read coverage. mpileup doesn't work with PAF file so I am looking to convert PAF format to SAM/BAM format. If there are alternatives to mpileup, I am open to the suggestions. Is this better ?

It's really weird that they provide sam2paf but not paf2sam. Could it be that the conversion in that direction isn't possible? Perhaps paf is missing information that would be essential for sam?

@terdon I am not a specialist but it's stated that CIGAR can take different format depending on the input/output options. It might also be not written in the PAF in some case. Is that CIGAR essential to sam ?

I would take Heng Li's advice.  minimap2 is fast. Just add the -a option when you re-run, which outputs the alignment in sam instead of paf

The paf format does not have the CIGAR string and base quality information that are typically present in a BAM/SAM file. While it is in principle possible to write a paf2sam converter, it's going to be very difficult to regenerate the CIGAR string with all the indel information (and maybe it's probably impossible to do so) that you'll need in a BAM/SAM file. If your motivation is to get coverage information, it might make more sense to convert the paf file into a bed file, before using a tool like bedtools to convert it into a coverage track.

note that PAF files can have a "cg:" tag that is a cigar string. it is an optional output from minimap2 for example (there is also cs: tag that is a different encoding of the sequence differences, but cg may be more common)