cannin/biostars_qa huggingface dataset:
Question: Successfuly install MaSuRCA

> tar xfz MaSuRCA-3.2.1.tar.gz
>
> cd MaSuRCA-3.2.1
>
> sudo ./install.sh

Configure and run config.txt

> MaSuRCA-3.2.1/bin/masurca MaSuRCA-3.2.1/bin/config.txt

Verifying PATHS...
jellyfish OK
Error: runCA not found or failed to run

Help me to get through this error.

Answer: Solved:

> PATH=~/tools/MaSuRCA-3.2.1_12062016/CA/Linux-i686/bin:$PATH

> PATH=$PATH: ~/tools/MaSuRCA-3.2.1_12062016/CA/Linux-i686/bin

> sudo env "PATH=$PATH" ~/tools/MaSuRCA-3.2.1_12062016/bin/masurca ~/tools/MaSuRCA-3.2.1_12062016/bin/config.txt
Question: Hello,

I have a tsv file with 1s and 0s. It has about 20 rows and 100 columns. I have created a dataframe. Now, I want to create a plot/graph which can show a particular colour/pattern for all 1s and another for the 0s with rownames in x-axis and column names along the y-axis. Is this possible in **R** ? I tried creating normal plots but this gives error saying "rows and columns are not equal".

Thanks in advance.

Answer: Simplest way, using *image* plot:

image(x, col = c("red", "blue"), main = "base image plot")

And for more flexibility, we can try *ggplot*:

library(ggplot2)

#convert wide to long format
plotDat <- reshape::melt(x)

#and plot
ggplot(plotDat, aes(X1, X2, fill = c("red", "blue")[ value + 1 ])) +
geom_tile() +
scale_fill_identity() +
theme_minimal() +
ggtitle("ggplot geom_tile")

Example data

set.seed(1); x <- matrix(sample(0:1, 25, replace = TRUE), nrow = 5)
colnames(x) <- paste0("col", 1:5)
rownames(x) <- paste0("row", 1:5)
x
# col1 col2 col3 col4 col5
# row1 0 1 0 0 1
# row2 0 1 0 1 0
# row3 1 1 1 1 1
# row4 1 1 0 0 0
# row5 0 0 1 1 0
Question: im newbie to NGS and linux. can anyone kindly let me know the easiest way to split the fasta file containing chromosome and scaffold into individual files i.e files containing one chromosome per file...
thanks

Answer: A modification of Pierre's answer from here : https://www.biostars.org/p/76716/

awk '/^>/ {F=substr($0, 2, length($0))".fasta"; print >F;next;} {print >> F;}' < ref_genome.fasta

Just replace `ref_genome.fasta` with your file
Question: Hi everyone,

For an allele mining project I need to generate in silico NGS data. I have chosen to use Grinder
to do so. I would like to add randomly generated SNPs to my reference genome, is any software
available to do that ?

Thanks for help!

Answer: You could try `mutate.sh` from `BBTools/BBMap`. `mutate.sh` can also produce a VCF file of the SNPs added. You can also generate random Illumina reads with the BBTool `randomreads.sh`
https://sourceforge.net/projects/bbmap/
Question: I'm trying to merge different bam files with the same name (and, theoretically, the same sequence). They are the same set of libraries sequenced at two different times, and I want to merge the matching pairs in order to simulate greater sequencing depth.

I can do this just fine with two bams using `picard MergeSameFiles` or `samtools merge`, but the issue is that I have 96 bams in each folder. I'd like to do this programatically, not manually, and be able to reproduce the process with different datasets in the future.

My hunch says that the simplest way to do this would be to use the bam filenames: loop over my two folders, find the bams that share a filename, and merge each pair into a single output file, but my shell chops are still rough and I am hitting a wall.

I've been starting with getting a list of bams of interest (I've also attempted dumping this list to a file with `-fprint0`):

find_bams() {
find "$run_folder" -type f -name "*.bam"
}

Then, what I think I should do is loop over the list:

for i in $(find_bams)
do
s=$(basename "$i" .bam)
picard MergeSamFiles I="$i" O="$s".bam
done

This is where I get stuck. First, there needs to be *two* input files to merge, and second, those two inputs must have matching filenames.

This is more of a shell scripting problem than a bioinformatics software problem, but I imagine I'm not the only one who has had to do this. Any help would be greatly appreciated.

**Edit**: Solved with the help of [h.mon][1] and [finswimmer][2]; see comments

[1]: https://www.biostars.org/u/6093/
[2]: https://www.biostars.org/u/37605/

Answer: Unless you want to a write a more general script to tackle the problem in a more robust fashion, I would try a simpler approach. If the sets of bams have exactly the same names (as you said they have), I would do something like:

cd /path/to/bams1
for i in *.bam
do
picard MergeSamFiles I=/path/to/bams1/"$i" I=/path/to/bams2/"$i" O=/path/to/merged/"$i"
done

If you keep the same structure of folders and filenames for future datasets, this simple solution will work just fine.
Question: Hi

I would like to plot RNASeq data that I have downloaded from TCGA in a PCA plot. I have found some great guides on how to plot the actual data in PCA using r in ggplot2 and such but my main question is what format data should I plot?

I currently have raw counts and RSEM data. Should I input raw counts into something like edgeR or deseq2 and filter for expression by cpm first? Should I normalise it? Should I stabilise variance using rlog2? Or convert to TPM and plot that? Argh I'm so confused. Grateful for any advice you can give me :)

Answer: I would recommend you to input your raw counts into DESeq2, run the pipeline, convert normalized reads to rlog (regularized log transformed reads) and then just run the plotPCA function from DESeq2. It is very easy if you are familiarized with the program.

# Creating deseq2 object
dds <- DESeqDataSetFromMatrix(countData = inputData,
colData = samples,
design = design)

dds <- DESeq(dds, betaPrior = betaPrior)
# Regularized log transformation for different analysis (clustering, heatmaps, etc)
rld <- rlogTransformation(dds)
pca <- plotPCA(rld, intgroup = c(colGroups))

The idea behind using rlog transformation for Quality Control checks is described in [DESeq2 paper][1]: "[...] Therefore, we use the shrinkage approach of DESeq2 to implement a regularized logarithm transformation (rlog), which behaves similarly to a log2 transformation for genes with high counts, while shrinking together the values for different samples for genes with low counts. It therefore avoids a commonly observed property of the standard logarithm transformation, the spreading apart of data for genes with low counts, where random noise is likely to dominate any biologically meaningful signal[...]"

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
Question: Hi,

I am trying to run PFAM locally. I need to download some databases from the website: ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/ But there are too many files and I am not so sure which to download. All I want to do is to get the multi-domain architecture of my fungal protein sequences. Does anyone have any idea which databases I would need to download?

Thanks in advance!

Answer: Assuming you want to run hmmer on it: ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz.
Question: Hello everyone!

I'm relatively new to Bioinformatics so please bare with me. A Ph.D student in my lab has asked me to analyze ChIP-Seq data and determine whether peaks fall into exons, introns, or exon-intron junction categories. The gene she is looking at is TNFAIP3 and the mark she wants me to analyze is Pol II. I am using the hg19 human model.

Unsure how to start I am thinking of doing the following:

1. Download BED files for TNFAIP3 and get individual files for Exons, and Introns.
2. Annotate peaks of BED files using a python script already written using Homer.

After this I am unsure how to continue. Am I done once I have these annotations in a txt format that can be opened on excel? Her second question revolves around analyzing peaks for common sequence motifs so how would I prepare this information to continue into that?

Answer: Since you have a NarrowPeak file of ChIP signal peaks, and you have aligned ChIP signal, you can use either file (with slight differences) to measure whether one feature (exon, intron, exon-intron boundaries, intron-exon boundaries, 5' vs 3' exons) tends to have a greater signal (normalized to feature length). Bedtools should be able to do this, once you've constructed your different feature files in BED format (or similar), and I'm sure there are more sophisticated solutions.
Question: Recently, I searching paper about time series or time course, I found that these two words equally used by researchers. However, time series seems used more wider. wikipedia also use time series (https://en.wikipedia.org/wiki/Time_series). What the difference between these two phrases? Is there any specific context that one phrases more frequently used?

Best,

Pengcheng

Answer: @Jean-Karim may be clear enough, but let me try giving more details. I got confused about these terms sometimes as well. Le's check dictionary first. The relevant definition are:

> **course**
/kɔːs/
*noun*
noun: **course**; plural noun: **courses**
**1**. the route or direction followed by a ship, aircraft, road, or river.

> **series**
/sɪəriːz,-rɪz/
*noun*
noun: **series**; plural noun: **series**; modifier noun: **series**
**1**. a number of events, objects, or people of a similar or related kind coming one after another.

So, "time course" means course of the time, like in "course of the road". With out an "experiment" the term "time course" does not have any data and therefore we can't analyze anything. That's why we say "***time course experiment***". And the results of a time course experiments to be meaningful, they must be series of data points. Otherwise, why do we care doing time course at all? So we get "***time course experiment data series***". That is too long, let's call "***time series data***" for short :-)

In other areas of data analysis, weather forecast or stock market for example, they have time series data too. However, the "experiment" is measurements.

With that said, and strictly speaking, "***time series data***" is not meaningful if you take it literally because there is no "series of time" alone. It must be understood as a shorthand of ***time course experiment data series***.
Question: I have a total of 300 million reads in my R1 and R2 raw fastq files.

How come fastQC count about 50% duplicates in each .fastq files, but after alignment with STAR (default parameters), none of the reads have the 1024 flag (PCR or optical duplicates) when doing `samtools view -c -f 1024 mydata.bam`?

Answer: Because STAR with default settings doesn't mark duplicates, you have to use `--bamRemoveDuplicatesType UniqueIdentical` to have duplicates marked (or `--bamRemoveDuplicatesType UniqueIdenticalNotMulti`).

**edit:** also, to elaborate on the point https://www.biostars.org/u/18713/ makes:

STAR (if you choose to mark duplicates with the correct parameter, which apparently you didn't) and FastQC use very different methods of detecting duplicates, so it shouldn't be surprising if the number os estimated duplicates is different. The [documentation for FastQC][1] states:

>To cut down on the memory requirements for this module only sequences which first appear in the first 100,000 sequences in each file are analysed

And also:

>Because the duplication detection requires an exact sequence match over the whole length of the sequence, any reads over 75bp in length are truncated to 50bp for the purposes of this analysis.

Whereas STAR uses mapping position:

> UniqueIdentical
>
> mark all multimappers, and duplicate unique mappers.
>
> The coordinates, FLAG, CIGAR must be identical

[1]: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/8%20Duplicate%20Sequences.html
Question: Dear biostars,

Is there a statistical method to tease whether an elevated polygenic risk score is independent to a monogenic risk? I have a cohort of patient and controls who have undergone whole genome sequencing. I have then done collapsing rare variant testing using SAIGE-GENE and found a gene that is strongly enriched in the disease cohort. I have then applied a validated polygenic risk score (PRS) and again found a statistically significant enrichment in cases over controls both including and excluding those cases that are part of the monogenic hit.

A subset analysis of the PRS cases and controls with variants that qualified for the monogenic hit show a near doubling of the PRS in cases over controls but the numbers are likely too small to reach statistical significance.

My question is whether there is a method to see if the PRS and monogenic hits are statistically independent from each other. My hypothesis is that those with monogenic risk also have an elevated PRS and it is that which is what pushes them into having the phenotype as the OR of having disease in the presence of the monogenic variants is 2.5 with a penetrance of 0.28.

Many thanks for your time

Answer: "Statistical independence" is really a heavy lift; and the only physical way to achieve this would be to have the gene independently **inherited** from all of the other polygenic risk variants; and the only way I can think of for this to happen would for the gene to be on its very own chromosome. Otherwise there will always be mutations in LD with your gene that capture some fraction of the risk. However the **correlation** could be low (if, say, most of the PRS is coming from other loci).

Have you tried simply plotting the burden score (x-axis) and PRS (y-axis)?
Question: Hi all,

Currently, I am working on my very first RNAseq study and have met a dilemma where inputs from more experienced bioinformaticians would be amazing.

For a differential gene expression study in a non-model organism, a de novo reference transcriptome was assembled from 300 M reads in Trinity. For 3 experimental conditions (1 negative control, 1 positive control and the treatment of interest) triplicate samples were sequenced with a depth of 25 M reads.

The reference transcriptome was annotated with Trinotate.

For differential gene expression determination the Kallisto/Sleuth pipeline was being used - and here comes my dilemma of best practices:

A number of the Trinity transcripts could not be annotated by Trinotate (NA) and is being dropped in the Sleuth analysis when using the *"so <- sleuth_prep(s2c, ~treat, target_mapping = annotation, aggregation_column = 'gene')"* expression.

I played around with the annotation file and replaced the NA's in the gene column with the corresponding Trinity transcript IDs, which included some of the transcripts as significantly differentially expressed.

**What is the right thing to do?**

- Would you let Sleuth drop the non-annotated transcripts, even though some of them are significantly differentially expressed?

- Or, would you include these transcripts in the "gene" column with their corresponding Trinity transcript IDs, even though they cannot be analyzed on the gene-level (the transcript isoforms cannot be collapsed in the analysis like the annotated ones)?

Thank you!

Cheers, Birgitte

Answer: This is a tricky one, and often "the right thing to do", becomes "what I can logically defend". Frequently in research, you'll reach a point where there's not really the right way to do something, but the most sensible for your current situation that you can defend when asked.

Specifically to your question, I think it's logical for you to change your annotation file so that you include the respective trinnotate IDs where there is an NA. Just because trinity can't annotate something that it thinks is a transcript, doesn't mean that it's not potentially interesting. Also try to keep in mind that you're trying to do an incredibly complex task, with a technology that is sub-optimal for its purpose (assuming short read sequencing), and this assembly will not be perfect. You're going to see noise, i.e. transcripts that trinity thinks are real when they're not (false positive), and vice versa (false negatives).

The most logical thing to do depends on your experimental question, but if it's "Whats different in the assembled transcriptome between my conditions?", then focus on those that are differentially expressed. Anything where you've got a useful trinnotate annotation is fine, but those that have your replaced trinnotate IDs might be worth a follow up. Consider extracting one of those sequences and blast searching to inspect the results in a bit more detail.
Question: Hello,
Following are the list of dual indexed illumina libraries we have used for illumina run.But we were not able demultiplex the data by using the barcodes.Because 99% of the data going as undetermined.

Index1 ------ Index2 ( case1 )

ATGAGC ---- TGAACCTT

ATTCCT ---- TGAACCTT

CAAAAG	---- TGAACCTT

CAACTA	---- TGAACCTT

CACCGG ---- TGAACCTT

There after we tried demultiplex with reverse compliment of index 2 and No changes in Index1(list are below ).This worked well and able to determine 98% of the data.The question is why the demultiplex only worked with reverse compliment of index2 and not in the case1(without reverse compliment)

Index1 ------ Index2 (reverse compliment)

ATGAGC ---- AAGGTTCA

ATTCCT ---- AAGGTTCA

CAAAAG ---- AAGGTTCA

CAACTA ---- AAGGTTCA

CACCGG ---- AAGGTTCA

Answer: The index2 (i5) sequence orientation is dependent on the sequencer used. Your original i5 sequence for barcode A501 is appropriate for a HiSeq 2000/2500, MiSeq, or NovaSeq, but would need to be reverse-complemented for MiniSeq, NextSeq, or HiSeq 3000/4000 (you can see this [in the following PDF from Illumina](https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/experiment-design/illumina-adapter-sequences-1000000002694-06.pdf)). This turns out to be rather annoying if you run a sequencing facility with multiple types of sequencers and people provide the barcodes they use. You then either have to check the orientation manually (that's annoying) or what for failures like this and reverse-complement accordingly (this actually works OK most of the time). I ended up writing a python package for [our demultiplexing pipeline](https://github.com/maxplanck-ie/TheWhoTheWhatTheHuh) that tries to determine the proper orientation of the barcodes in each lane [by looking at what's in the BCL files produced by the sequencer](https://github.com/maxplanck-ie/TheWhoTheWhatTheHuh/blob/master/pyBarcodes.c#L18-L40).

Update: Please see the comment below from https://www.biostars.org/u/3616/ for an update on the NovaSeq!
Question: What is the appropriate way for determining that the application and workflow licenses for NGS apps can be used (without any IP infringement) to host the apps/workflows on a pay-per-use infrastructure server in a localized region?

Answer: Much the same as any software. Each tool will be available under a particular license. If they're open source that typically means you are free to make money on the software (though there are some licenses which do not allow this IIRC). For tools which are not free/open-source, there will often be 2 license options, a commercial one, which enables you to charge for its use subsequently, and an individual/research/acadmic type license which will typically be cheaper but you cannot use 'for hire or reward'.

IP is a somewhat different matter, as it depends on the specific nature of what you're building from it as to whether it would infringe IP.

You will likely need to give us more information about where you are/plan to use/deploy the software, and what you plan to do, but long story short you need to consult the licenses for each bit of software you plan to use.

*The usual I am not a lawyer disclaimer applies.
Question: I have a fasta file like below, I like to seprate only coordinate and add one column with "hello" name.

>Human|chr16:80372593-80373755 | element 4 | positive | neural tube[6/10] | hindbrain (rhombencephalon)[10/10] | midbrain (mesencephalon)[10/10]
gtgaCAGAGACAGACAGTGACAGAGACAgattttagaatttgaacaaaggtaaataagag

>Human|chr16:78510608-78511944 | element 12 | positive | hindbrain (rhombencephalon)[9/11] | forebrain[9/11]
AAGCTAGCTAATTGCTTCTTCAGTTGaagacctaaatgagttttaaagtgaaatgcatat

Expect file:

chr16:80372593-80373755 hello
chr16:78510608-78511944 hello

Answer: awk -F '|' '/^>/ { print $2 "\thello"}' jeter.fa
Question: I want to create two random subsamples from a BAM file. I don't want to just do "samtools view -s 0.5" twice because I want to put half of the reads into one file and half in the other, without replacement. Is there an easy way to do this? I've considered using samtools to output the alignments in SAM format and pipe to the shuf and split commands, I think I'd need to do some extra things to keep the paired-end reads together. Anyone have any ideas?

Answer: Assuming single-end reads, something like this should work:

samtools view -H foo.bam > f1.sam
cp f1.sam f2.sam
samtools view foo.bam | awk '{if(NR%2){print >> "f1.sam"} else {print >> "f2.sam"}}'

I'm sure there's a one line version of that, but this is simple enough.

For paired-end reads, just change `NR%2` to `NR%4<2`. Again, there's likely a bug in there, but it's something to start with.
Question: I am planning to merge two FastQ file from NGS.

And I have question with FLASH algorithm.

I have below File Forward and reverse

1\. Forward

```
@CP000143_994500_994663_0:0:0_0:0:0_0/1
GCTTCTGCGACCGCGCCCTCGTCGTCTACCGCGGCACGCTGAACGGCGAGTTCGCGGGCGAGACGCTCGACAGCGACCTGCTCCTGGCCGCCGCCTCGGG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
```

2\. Reverse

```
@CP000143_994500_994663_0:0:0_0:0:0_0/2
GCTCTCGCTGGAGGACGGGGACAGGGCCATGGTCATTCAGGCCTCCTTTCGTTGGGCCCGCGCGCCCGAGGCGGCGGCTAGGAGCAGGTCGCTGTCGCGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
```

If Reverse Complement and Merge together results will be.

GCTTCTGCGACCGCGCCCTCGTCGTCTACCGCGGCACGCTGAACGGCGAGTTCGCGGGCGAGACTCGCGACAGCGACCTGCTCCTAGCCGCCGCCTCGGGCGCGCGGGCCCAACGAAAGGAGGCCTGAATGACCATGGCCCTGTCCCCGTCCTCCAGCGAGAGC

Compare both Sequence

```
GCTTCTGCGACCGCGCCCTCGTCGTCTACCGCGGCACGCTGAACGGCGAGTTCGCGGGCGAGACGCTCGACAGCGACCTGCTCCTGGCCGCCGCCTCGGG
^----------------------------------^
TCGCGACAGCGACCTGCTCCTAGCCGCCGCCTCGGGCGCGCGGGCCCAACGAAAGGAGGCCTGAATGACCATGGCCCTGTCCCCGTCCTCCAGCGAGAGC
^----------------------------------^
```

And There is 3 mismatch during merge two sequences

But that mismatch part Results will follow reverse complement. Anyone has idea why it happen like this?

Thank you!

Answer: In the case of mismatches between the two sequences, FLASH will use the nucleotide with the highest base quality. Of course with fake qualities you're going to get weird results on occasion. I should note that mismatches like this only happen is something like 0.01% of alignments in bisulfite sequencing datasets...so in RNAseq it should be even less common (so it shouldn't be common enough to care about).

BTW, in the source code for FLASH (combine_reads.c, around line 400), this is how mismatches with equal phred scores are handled:

```
if (*qual_1 > *qual_2) {
*combined_seq = *seq_1;
} else if (*qual_1 < *qual_2) {
*combined_seq = *seq_2;
} else {
/* Same quality value; take the base from the
* first read if the base from the second read
* is an 'N'; otherwise take the base from the
* second read. */
if (*seq_2 == 'N')
*combined_seq = *seq_1;
else
*combined_seq = *seq_2;
}
```

One could argue about their choice there at the end, but that's why you get the results that you do.
Question: Hi All,
I intersected to data sets and I have two columns like:

m01 m01
m01 m02
m01 m05
m01 m032
m01 m02
m01 m01
m02 m06
m02 m01
m02 m02
m02 m09
... ...
m0500 m023
...

I would like to get number of matches with others like:

m01 m02 m03 ......

m01 25 45 98 .....
m02 90 223 12 ......
.
.

Would you please help me how I can do that?

Thank you very much

Answer: You won't get a matrix, but the answer would be the same.

sort 2-col.txt | uniq -c

Or with awk.

awk -v OFS="\t" '{cnt[$0]+=1} END{for (key in cnt) {print key,cnt[key]} }' 2-col.txt

Or with perl.

perl -lane '$cnt{$_}++; END { foreach $key (keys %cnt) {print "$key\t$cnt{$key}"} }' 2-col.txt