cannin/biostars_qa huggingface dataset:
Question: Hi All,

I am wondering what is the most widely used method/program/tool to automate the launch of bcl2fastq program as soon as the run finishes off on the Illumina machine. I am about to start writing my own custom shell script for the same; but before I do that I want to know whether ready-made solutions are already available.

The run folder nomenclature is as follows:

**My idea would be like this -**

Keep looking for folder containing the "Instrument Serial No." in its name:
for every 15 minutes :
Look inside the run log file for something which says - sequencing complete
if YES:
Launch bcl2fastq
else PASS:

If someone has better idea which could be implemented, please do let me know.

Answer: Hi

Coming back to my own question with an answer. After a lot of careful considerations and talking to the Illumina people, it turned out that generation of `CopyComplete.txt` (NovaSeq) is a good trigger to start bcl2fastq. Here is a little bit more information on that:

1. `SequenceComplete.txt` indicates that the sequencing run has finished.
2. `CopyComplete.txt` is created by the Universal Copy Service (UCS) when all files have been copied to their destinations and run completion signal has been triggered.
3. `The RTAComplete.txt` file indicates that the images that are generated by the system have been converted to Basecalls. The Basecalls are stored as .bcl files. These bcl files can then be used as input for BCL2Fastq to produce the fastq files.
Question: I have a new strain of Bacteria for which no reference genome is available at NCBI.
The strain is sequenced using Illumina and the paired-end reads were obtained.

The pair-end reads were taken and subjected to De-novo assembly using CLC genomics workbench. Once the process of de-novo assembly was complete, the following results were obtained.

Name	Consensus length	Total read count	Single reads	Reads in pairs	Average coverage
contig_1_mapping	30575	24185	749	23436	116.3326574
contig_2_mapping	4087	2703	111	2592	97.23562515
contig_3_mapping	5764	3835	187	3648	97.76040944
contig_4_mapping	7302	3488	152	3336	69.38989318
contig_5_mapping	26530	16489	553	15936	91.15295891
contig_6_mapping	5576	3883	119	3764	102.2159254
contig_7_mapping	42110	33820	1166	32654	117.8527666

Then contig measurements were also obtained like `N75	N50	N25	Minimum	Maximum	Average	Count	Total`

And the following contigs details

Name	Modified	Description	Size	Linear
contig_1	Mon May 15:35:29 IST 2017	Average coverage: 116.33	30575	Linear
contig_2	Mon May 15:35:29 IST 2017	Average coverage: 97.24	4087	Linear
contig_3	Mon May 15:35:29 IST 2017	Average coverage: 97.76	5764	Linear
contig_4	Mon May 15:35:29 IST 2017	Average coverage: 69.39	7302	Linear
contig_5	Mon May 15:35:29 IST 2017	Average coverage: 91.15	26530	Linear
contig_6	Mon May 15:35:29 IST 2017	Average coverage: 102.22	5576	Linear

Now how do I analyze and annotate these results which were obtained after doing de-novo assembly of the strain? Is there any NCBI tools to annotate these contigs?

Answer: Initial annotation:

- prokka

- [RAST][1]

Then you may wanto to use antiSMASH.

You may want to ckeck for contaminations with blobtools.

Read [this][2] paper.

[1]: http://rast.nmpdr.org/
[2]: https://microbialinformaticsj.biomedcentral.com/articles/10.1186/2042-5783-3-2
Question: The following is the output of samtools flagstat command on bam file (paired-end) generated after markDuplicate of Picards.

7417232 + 0 in total (QC-passed reads + QC-failed reads)
287618 + 0 duplicates
4534962 + 0 mapped (61.14%:-nan%)
7417232 + 0 paired in sequencing
3708616 + 0 read1
3708616 + 0 read2
4528278 + 0 properly paired (61.05%:-nan%)
4534962 + 0 with itself and mate mapped

I am having difficulty in understanding whether the duplicates are pairs or single. If there are total of 7417232 pairs and out of them 287618 pairs are duplicates means, there are 3% of duplicate reads in my data. Is my understanding is correct ?

Answer: There aren't 7417232 pairs, but total reads. There are 3708616 actual pairs of reads. It's easiest to figure out exactly what's going on by just looking at the C code in bam_stat.c (it's a short file). The 287618 + 0 duplicates numbers are incremented every time a read is marked as a duplicate. So, if there are 100 duplicates for one read then that number will increase by 100 (not 1). Practically speaking, that would mean that there are 143809 duplicate pairs, which could all be the same duplicate (or not). Regardless, that's about 3% (closer to 4% I think) of duplicates, so that part of the interpretation is spot on.
Question: I have single-end fastq.gz files.

when I run RSEM with `--gzipped-read-file` it says invalid number of cammands and when i dont use it, it says `FATAL ERROR, exiting. What should I do?

rsem-calculate-expression --star -star-path /share/pkg/star/2.5.3a -p 4 EH_25_data/EH025x045_L001/L001_R1_001.fastq.gz src_eh/rsem-ref/ref ./EH_25

Answer: Use --star-gzipped-read-file instead of --gzipped-read-file !
Question: I would like to verify if all the protein in specific pathways (e.g. secretion systems, two component regulatory system) are present in a genome.

Is there a way to do this automatically iterating the search for many genomes?

Answer: https://gist.github.com/RamRS/d772a4e8e562a45f0207e3ea1d18472e

Cleaned up the post and moved to a Gist by https://www.biostars.org/u/8494/ on 22-Apr-2022.
Question: Hello,

I have a file -

```
GeneID sample1 sample2 sample3 sample4 gene_length
gene1 count1 count2 count3 count4 length1
gene2 count1 count2 count3 count4 length2
gene3 count1 count2 count3 count4 length3
gene4 count1 count2 count3 count4 length4
```

When imported to R using `read.table()` I get columns assigned by default as V1, V2, V3 etc, but how can I set the first row as the column names?

GeneID sample1 sample2 sample3 sample4 gene_length

Thanks

Answer: Add Argument `header=TRUE`
Question: Hi there,

I'm sorry for the newbie question but I'm confusing and getting lost in my thoughts.
Having `WES` data, until the moment I have been using `DNASeq` aligner such as `BWA MEM`. But I'm wondering whether it makes sense to use `RNASeq` mapping algorithm (`STAR`...) or not and why.

For other hand lets say that I want to consider the alignment of the reads against only the exons of some genes. How this would affect the aligner algorithm decision?

Thank you in advance,

Answer: Stick with a DNA based aligner like BWA (recommended in the [GATK best practises][1]). DNA sequencing (such as WES), is sequencing DNA fragments pulled down by a pre-designed capture kit, whereas RNA sequencing is sequencing the actual mRNA. When aligning RNA reads, the aligner (STAR, HISAT2, etc) will expect reads to split across exons, as such the appropriate penalty for that is tuned - that sort of event shouldn't happen with DNA.

[1]: https://software.broadinstitute.org/gatk/best-practices/bp_3step.php?case=GermShortWGS
Question: Hallo,

I have some problem with importing data using GEOquery. I download an ExpressionSet for some further analysis, however in my record only few GSMs from all are interesting for me. My question is how to modify my code to download only the selected GSMs and analyze it further with limma package. I am grateful for any suggestions.

eset<-getGEO('GSE13937')[[1]] #-----> but I would like the data only for e.g. GSM26056, GSM26076, GSM26102, and GSM26123

And, then perform a standard analysis like:

```
design <-cbind(normal=c(1,1,0,0),
tumour=c(0,0,1,1))
c <-makeContrasts(normal-tumour,levels=design)
fit <-lmFit(eset, design)
fit2 <-contrasts.fit(fit,c)
fit2 <-eBayes(fit2)
glist <-topTable(fit2, n=50, adjust="BH", sort.by="P", p.value=0.01)​
```

Answer: Well, subsetting an ExpressionSet is pretty easy:

eset1=eset[,c("GSM426056", "GSM426076", "GSM426102", "GSM426123")]

However, you might not want to do things this way. Instead, you can create your design matrix with all "0" except for the two "normals" and the two "tumors". The rest of your code can remain unchanged (the contrasts, testing, etc.). You get the benefit of a more robust measure of variance if you include all samples even though the testing will be based on only your four samples.
Question: Anybody could help me explain values below? What I know is they are from the results of edgeR package and used for logFC calculation. I am a fresh man in R and thank you very much!

```
id logFC logCPM PValue FDR A1 A2 A3 B1 B2 B3
c82096_g4_i5 -9.07576 5.512818 2.87E-54 4.46E-49 50.32 38.61 29.09 0 0 0.18
c76837_g1_i2 -12.4923 4.185135 5.83E-50 4.53E-45 12.22 11.85 11.96 0 0 0
c74136_g1_i4 -12.7555 4.456949 3.30E-49 1.71E-44 29.37 18.52 16.64 0 0 0
c82993_g3_i5 12.43865 4.169207 5.75E-47 2.23E-42 0 0 0 3.09 2.95 2.46
```

Answer: The updated numbers make more sense. In this case, fold-changes like those you posted look reasonable. I'm more familiar with DESeq2, but in general these numbers are the result of maximum likelihood expectation after shrinking the variance. Once you start approaching infinity, the standard error of those values start sky rocketing, so you can take any values greater than 6 or 7 (i.e., a fold change greater than 100 or more) as "near infinite", at least unless they have some meaningful (i.e., non-zero) counts across groups.

BTW, you'd likely get smaller fold-changes with DESeq2, since the large values you're seeing are partly due to having only 3 samples per-group (i.e., they're not that reliable and would, thus, get shrunken a good bit, though they'd still be quite large).

Edit: I probably should have written "profile likelihood" above instead of "maximum likelihood", since I expect dispersion is just factored out.
Question: I need a way to assign a number to genes representing how broken it is depending on how many SNPs with minor alleles appear on the gene. Based on just a handful of SNPs though, not every SNP.

So lets say for example I know a few SNPs for COMT:

```
Rs165599 AA
Rs6269 AA
Rs737865 AA
Rs4633 TT
Rs769224 GG
Rs4680 AA
```

and some SNPs for STAT3:

```
Rs744166 GG
Rs8069645 AG
Rs9891119 CC
Rs6503691 CC
```

and I wanna find out which gene is more defective (as in which has the most diminshed ability to encode the corresponding protein), how can I do that? These genes obviously have way more SNPs than just those ones, and in every gene I will have many minor alleles so I'm guessing its only certain SNPs that cause problems for the gene. So could I find a number that tells me how much a bad allele of a particular SNP breaks the gene and then use that number to get an overall idea of how malfunctioning the gene is?

Answer: In addition to what has already been suggested there is also CADD and RVIS that you should check out, CADD is much like SIFT and PolyPhen, RVIS is a statistic calculated that estimates how tolerant/intolerant genes are to mutations, which is probably very important to take into consideration for what you are looking in to.

The other thing to keep in mind is that not all mutations, even rare ones, will necessarily be detrimental to the protein's function. Some will have minor impacts, many will be essentially neutral in effect, and some may even enhance the protein's function. These are all important concepts from molecular evolution that aren't taken into account as often as they should be by many human geneticists.
Question: Hi,

I have a list of Ensembl protein pairs(homologs) with their respective ensembl protein IDs. I want to find a way to convert these Ensembl protein IDs to Uniprot IDs. for example given: [ENSP00000361930][1] I want to get: `1433B_HUMAN`

Is there a function in biopython that does this conversion? I really need to script it because the list is huge and I cant do it manually.

Thank you

[1]: http://www.ensembl.org/id/ENSP00000361930

Answer: use uniprot ID mapping http://www.uniprot.org/uploadlists/
Question: Dear all,

in the light of a new article in PLOS : "Evaluation of Nine Somatic Variant Callers for Detection of Somatic Mutations in Exome and Targeted Deep Sequencing Data" (http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0151664) that recommends the use of Mutect, Strelka and Virmid for somatic mutation calling, wondering if you guys can share your experience with these 3 algorithms, or if you have any other recommendations. Thanks a lot,

bogdan

Answer: I always interpret some of these results with a healthy grain of salt, in particular because the actual truth set of mutations usually isn't defined, so that there is no good analysis of Precision, Recall, and other measures. The authors also generally assumed that the four callers that had tended to converge on one another, and had the most agreement with each other, where also the correct callers. I happen to agree with them, but that kind of reasoning can be flawed if they are suffering from systematic bias due to a similar calling method. That said, of the recommended callers I've only used MuTect extensively. I'm somewhat surprised that FreeBayes wasn't included as it is widely considered to be an excellent Haplotype based caller as well. MuTect, in my experience, works quite well, although I've found getting the old version is now somewhat of a pain to find and download, and get good documentation for. There are also issues I believe with making sure you are using the correct Java version, which is a perennial problem for Java software in general and anything that uses the GATK in any way in particular.

I'm now looking at testing EBCall and Virmid in my own pipeline as they weren't in there. I've been using MuTect, FreeBayes, Scalpel, VarDict, Pindel, and Platypus.
Question: Hi,

I have obtained BLAST output file in TSV file format which looks like -

```
# BLASTN 2.2.29+
# Query: SI2.2.0_06267 Si_gnF.scaffold02592[1282609..1284114].pep_2
# Database: ./nucleotide/Sinvicta2-2-3.cdna.subset.fasta
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 1 hits found
SI2.2.0_06267 SI2.2.0_06267 100.00 240 0 0 1 240 1 240 1e-126 444
```

This output file has no visual alignment information (which is shown normally by aligning the query and hit sequences on top and drawing |'s in between to show matching). Is it possible to re-generate or create that from the information given by TSV output? If yes, what would be an efficient way to do that?

Please note that, I wish to use TSV output in first place.

Thanks,

Answer: It is not possible to recreate the alignment without recomputing it. This is because there is no information about the location of mismatches and gaps. If you need multiple output formats it is better to use ASN.1 binary output format and reformat to the desired output format with `blast_formatter`.
Question: Hi all,

I am having a hard time to extract the corresponding subgroup members in the heatmap generated by heatmap.2 in gplots.
I have looked here: http://stackoverflow.com/questions/15388082/understanding-heatmap-dendogram-clustering-in-r
and http://stackoverflow.com/questions/11489696/how-to-know-about-group-information-in-cluster-analysis-hierarchical
and here https://stat.ethz.ch/pipermail/bioconductor/2007-December/020521.html

I have a matrix m, the heatmap can be generated by

heatmap.2(m, colv=FALSE, dendrogram= "row", trace="none")

it gave me very good looking heatmap, now I want to extract the subgroup members.

I know I can use hclust and cutree to extract the information:

hc.rows<- hclust(dist(m))

plot(hc.rows) # it gives me just the dendrogram, and it looks the same as the dendrogram appeared in the heatmap, only differ in the order. it seems to me the order is reversed.

I can cut the tree based on either the number of group (k), or the height (h)

ct<- cutree(hc.rows, h=10) # it gives me 6 groups

rect.hclust(hc.rows, h=10) # draw red rectangles to mark the subgroups

table(ct)

ct

1 2 3 4 5 6

196 248 294 119 42 4

I can get the labels of each row, and the assigned the group number by:

tableclust<- data.frame(m,ct)

my question is that how do I map the cluster number here back to the dendrogram?
cluster 6 has 4 members, but it looks like to me cluster 6 appears in the most left of the dendrogram, cluster 5 has 42 members, it appears in the second most left (based on the width of the red rectangle in the dendrogram), cluster 4 has 119 members, but it appears not next to cluster 5....

I am confused.....

Thanks for your advice.

Tommy Tang

Answer: use `svg(width=50)` option. Increase it for bigger width. https://www.biostars.org/p/49786/
Question: When people ask for ld blocks on this forum, this link is often posted:

ftp://ftp.ncbi.nlm.nih.gov/hapmap/ld_data/latest

However, those files merely contain the correlations between snps. How do I use these to compute ld-blocks?

Would be nice if you could post an example script or invocation of a tool, not merely post a link :)

Thank you!

Answer: Method: https://academic.oup.com/bioinformatics/article/32/2/283/1743626

Data: https://bitbucket.org/nygcresearch/ldetect-data/src/master/EUR/fourier_ls-all.bed
Question: Hi there.

I have two bam files which followed the preprocessing including GATK.

and I would like to check these two bam files read quality as graph in order to understand easily.

Is any library or program that provides the read quality of these bam files containing with some charts?

Answer: RSeqQC might help you to get mapping statistics

bam_stat.py -i myFile.bam

for stats in graphical way

qualimap is nice one, you can find it [here][1] PS: I have never used it.

for other information this site would be useful for you

http://barcwiki.wi.mit.edu/wiki/SOPs/miningSAMBAM

[1]: http://qualimap.bioinfo.cipf.es/
Question: I want to create a .gmt file for visualizing Gene enrichment sets in the EnrichmentMap plugin for cytoscape

My input file is like this:

```
scigt000016 GO:0005515
scigt000021 GO:0005515
scigt000027 GO:0044464
scigt000010 GO:0005515
scigt000011 GO:0015074
```

And I want to convert it like this:

```
GO:0005515 NA scigt000016 scigt000021 scigt000010
GO:0044464 NA scigt000027
GO:0015074 NA scigt000011
```

So basically putting the GO-term in column 1, some random text in col2 and the genes from col1 in the input file on a line after the go-term.

I was thinking to use a for loop that greps the go-term from column 2 in the input file, then for each line append col1 to the end of the line. But really I am quite lost here.

Answer: sort on 2nd column, use awk to retrieve the 'previous' value of the 2nd column. If it isn't the same, start a new line.

```
sort -k2,2 input.txt | awk '{if($2!=prev) {printf("%s%s\tNA\t%s",(NR==1?"":"\n"),$2,$1);prev=$2;next;} else printf(" %s",$1);} END{printf("\n");}'

GO:0005515 NA scigt000010 scigt000016 scigt000021
GO:0015074 NA scigt000011
GO:0044464 NA scigt000027
```
Question: I have fasta file namely `119XCA.fasta` as shown below,

>cellulase
ATGCTA
>gyrase
TGATGCT
>16s
TAGTATG
I need to remove all the fasta headers, keep the sequences one by one and need to write file name as a fasta header. The expected outcome is shown below,

>119XCA
ATGCTA
TGATGCT
TAGTATG
I have used the following script `sed '/^>/d' foo.fa > out.fa` which remove the fasta headers but, i do not know how to manage to write file name as a header. Therefore, please help me to do the same.

Answer: Not the prettiest code in the world, but this will work.

Run it like so: `bash scriptname.sh /path/to/files/*.fasta`

for file in $1 ; do
cat $file | sed -e '1!{/^>.*/d;}' | \
sed ':a;N;$!ba;s/\n//2g' | \
sed '1!s/.\{80\}/&\n/g' | \
sed "s|>.*$|>${file##*/}|g" > $(basename "${file##*/}" ".fasta" ).fa
done

You can also do it as a oneliner for a single file if needed:

cat filename.fasta | sed -e '1!{/^>.*/d;}' | sed ':a;N;$!ba;s/\n//2g' | sed '1!s/.\{80\}/&\n/g' | sed "s|>.*$|>${file##*/}|g" > $(basename "${file##*/}" ".fasta" ).fa
Question: Hello everyone,

I am making a Ion Torrent pipeline by using tools from latest torrent suite. I am having hard time finding any sample data (with reference it was aligned against) using which I can do the tests (mapping, indexing and variant calling).

Are you familiar with any FTPs or something similar where I could find these?

Cheers!

Answer: There are two data sets from the NIST Genome in a Bottle project that would be good tests for Ion Torrent pipelines:

Ashkenazi trio: ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/AshkenazimTrio/

There is an ion_exome directory inside each of the trio members directories with Ion Torrent data.

NA12878 cell line: ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/ion_exome/

These are particularly nice because there is plenty of reference material on these samples so that you can validate your results.
Question: I use bcftools norm to split and left-normalize variants (like the line in file). After norm the variant is split correctly and annotated, my question is can the AF tag be used to filter out AF= below 20% in the command or is there another? So in this case the first split vatiant is filtered out. Thank you :).

before bcftools norm
chr4	187004216	.	TC	T,TT

after bcftools norm
chr4	187004217	187004217	C	-
AF=0.0305344;AO=130
chr4	187004217	187004217	C	T
AF=0.969466;AO=673

desired output
chr4	187004217	187004217	C	T
AF=0.969466;AO=673

Answer: You are certainly correct to normalise the file prior to filtering; otherwise, you would obtain unexpected results.

To filter on AF, use `bcftools view`, specifically the following parameters:

> -q/Q, --min-af/--max-af

Kevin
Question: I have a bunch of bulk mRNA sequencing pulled off of the TCGA. Feature names appear to be Ensembl gene IDs with a suffix. Here is an example:

> [995] "ENSG00000236246.1" "ENSG00000281088.1"
[997] "ENSG00000254526.1" "ENSG00000223575.2"
[999] "ENSG00000201444.1" "ENSG00000232573.1"

I am taking the intersection between these features and a set of Entrez Gene IDs. In order to do this I am using the biomaRt package to generate a mapping between Ensembl gene IDs and Entrez gene IDs. However, the only Entrez gene IDs I can find lack the suffixes. Here is the head of the table that maps Entrez genes to Ensemble genes:

entrezgene ensembl_gene_id
1 90529 ENSG00000001460
2 9235 ENSG00000008517
3 10747 ENSG00000009724
4 654364 ENSG00000011052
5 112611 ENSG00000013392
6 57210 ENSG00000022567

Can someone explain what the Ensembl suffixes mean and how to convert these names to Entrez? If this can be done with biomaRt, it would be ideal. Thanks.

Answer: The numbers are version numbers. There is information about stable ID versioning [here](http://www.ensembl.org/info/genome/stable_ids/index.html). You can just strip off the version numbers to use with biomaRt.