cannin/biostars_qa huggingface dataset:
Question: I have the following data:

Chromosome Start End XS XE TranscriptID GeneID Strand
36945 chr19 54754594 54756329 b'54754594,54756005,54756286,' b'54755063,54756202,54756329,' NR_110738 LOC101928804 -
36948 chr19 54769421 54771064 b'54769421,54770669,54771021,' b'54769891,54770937,54771064,' NR_110737 LOC101928804 -
36949 chr19 54769421 54771064 b'54769421,54770740,54771021,' b'54769891,54770937,54771064,' NR_110738 LOC101928804 -
36951 chr19 54785868 54835292 b'54785868,54816103,54834899,54835251,' b'54785899,54816541,54835167,54835292,' NR_110737 LOC101928804 -
36952 chr19 54785868 54835292 b'54785868,54816103,54834970,54835251,' b'54785899,54816541,54835167,54835292,' NR_110738 LOC101928804 -

Here you see transcripts from the same gene. Start is txStart and End is txEnd. For this gene (LOC101928804), can I say that the length of the gene is from 54754594 to 54835292 (from the start of first transcript to the end of last transcript). Or is this an oversimplification in some way?

Answer: In most all cases what you are describing should work to find the start/end of the gene. The way we build the refGene track is to BLAT all the NCBI sequences ourselves, though we also offer a track with all the NCBI provided alignmend: ncbiRefSeq. All of the NR_* accession numbers are experimentally validated sequences. the txStart is equal to the first position of exonStarts and txEnd equals the final exonEnds.

You can query these from our public MySQL server:

$mysql -h genome-mysql.soe.ucsc.edu -ugenome -Ee "select * from refGene where name2 like 'LOC101928804' and chrom not like '%alt'" hg38

*************************** 1. row ***************************
bin: 1002
name: NR_110738
chrom: chr19
strand: -
txStart: 54754594
txEnd: 54756329
cdsStart: 54756329
cdsEnd: 54756329
exonCount: 3
exonStarts: 54754594,54756005,54756286,
exonEnds: 54755063,54756202,54756329,
score: 0
name2: LOC101928804
cdsStartStat: unk
cdsEndStat: unk
exonFrames: -1,-1,-1,

And also pull out just the txStart and txEnd which may save you some time:

$mysql -h genome-mysql.soe.ucsc.edu -ugenome -e "select name,txStart,txEnd from refGene where name2 like 'LOC101928804' and chrom not like '%alt'" hg38

+-----------+----------+----------+
| name | txStart | txEnd |
+-----------+----------+----------+
| NR_110738 | 54754594 | 54756329 |
| NR_110737 | 54769421 | 54771064 |
| NR_110738 | 54769421 | 54771064 |
| NR_110737 | 54785868 | 54835292 |
| NR_110738 | 54785868 | 54835292 |
+-----------+----------+----------+

You can also pull out the selected info via point and click using http://genome.ucsc.edu/cgi-bin/hgTables.

Lou
UCSC GB
Question: I have a very basic question. I have created a local (mysql) database in my system. I would like to know what would be the best option for a front end. I had used html before and want to try out something new. I might also add some perl code in between for some computational work. Kindly help me in this regard..

Answer: Already Nice Answers: Here's my two cents.

I would recommend you to start with a basic template and add various components as required.

You can do the front-end -> database interaction in a simple way using [CGI.pm][1] with [DBI::mysql / DBD::mysql][2]

But the current state-of-art is to use [MVC][3] based frameworks (Catalyst is Perl based and Django is Python based MVC frameworks) for front-end with database interaction. A detailed list of Perl based frameworks are available [here][4]. Remember all of them are not capable of full MVC. For example [Template::Toolkit][5] will only take care of the view part. The latest in the block is [Mojolicious][6]. I am currently experimenting with this, got good reviews from fellow Perl fans. So if you want to try something different within Perl, you could give it a try.

To get the best look and user interaction for the interface you may need a bit of [CSS][7], [Java Script][8] and [Ajax][9]. If you are new to this, I would recommend you take a look at [W3C tutorials][10] and finally if you have enough time, you could learn a great deal of all these technologies from the online course videos : [Building Dynamic, Scalable websites][11] by [David J Malan][12].

[1]: http://search.cpan.org/dist/CGI.pm/lib/CGI.pm
[2]: http://dev.mysql.com/downloads/dbi.html
[3]: http://en.wikipedia.org/wiki/Model%E2%80%93view%E2%80%93controller
[4]: http://perl.apache.org/products/app-server.html
[5]: http://template-toolkit.org/
[6]: http://mojolicious.org/
[7]: http://en.wikipedia.org/wiki/Cascading_Style_Sheets
[8]: http://en.wikipedia.org/wiki/JavaScript
[9]: http://en.wikipedia.org/wiki/Ajax_%28programming%29
[10]: http://www.w3schools.com/
[11]: http://academicearth.org/lectures/building-dynamic-websites-http
[12]: http://www.cs.harvard.edu/malan/
Question: Hello!

Today I am working on concatenating paired files. In a directory, I have 300 paired files. Here's an example of the file names:
```
flL1_5495_L1PA6_reactivating_recurring_12_2bit_c.bed.fa.revcom.fa
flL1_5495_L1PA6_reactivating_recurring_12_2bit_plus.bed.fa
flL1_5495_L1PA6_reactivating_recurring_13_2bit_c.bed.fa.revcom.fa
flL1_5495_L1PA6_reactivating_recurring_13_2bit_plus.bed.fa
flL1_5495_L1PA8A_reactivating_recurring_03_2bit_c.bed.fa.revcom.fa
flL1_5495_L1PA8A_reactivating_recurring_03_2bit_plus.bed.fa
flL1_5495_L1PA8A_reactivating_recurring_04_2bit_c.bed.fa.revcom.fa
flL1_5495_L1PA8A_reactivating_recurring_04_2bit_plus.bed.fa
```

The file names are identical except for a family name (L1PA6 and L1PA8A in this example, but there are a few more), a level (12, 13, 03, and 04 in this example, levels range from 03 to 38 ultimately), and whether they are "plus" or "revcom". There is a matching revcom file for each plus file which I would like to concatenate. I have been working with a nested bash loop like so:

```
#!/bin/bash

subfamilies=( \
L1HS L1PA2 L1PA3 L1PA4 L1PA5 \
L1PA6 L1PA7 L1PA8 L1PA8A L1PA10 )

recurrence=( \
03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 \
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 \
35 36 37 38 )

for subfam in ${subfamilies[@]}; do
for level in ${recurrence[@]}; do
cat *${subfamilies[@]}_reactivating_recurring_${recurrence[@]}*.fa \
> ${subfam}${recurrence}.fa
done
done
```

I am trying to concatenate based on shared family names and levels, but my output is a mess. I get a lot of empty files. I'd like the output name to be simpler, something like "L1PA6_03.fa"

Maybe there's a better way to do this?

Thanks in advance!

Answer: How about something like:

```{sh}
for baseName in $(ls f* | sed 's/_2bit.*//' | uniq)
do
cat ${baseName}* > $baseName.fa
done
```

The idea is that the common part in the name of the files that have to be paired is everything before `_2bit`. So we list all names, truncate them with `sed`, and remove redundant entries. Enclosing this command in `$()`, we pass it as a plain space-separated string to `for`, and assign a `baseName` variable to loop on. The expression `${baseName}*` expand to the name of the two files to join together, and we reuse the `$baseName` to construct the output file name.
Question: I have this txt file by pasting different columns now how I can convert that to bed please because I need a bed then sorting that

I used

```
[izadi@lbox161 merged_asmribo]$ cut -f 1,2,3,4,5,6 main.txt > main.bed
cut: fields and positions are numbered from 1
Try 'cut --help' for more information.

FLP1 252 1523
RAF1 3271 3816
REP1 1887 3008
REP2 5308 6198
YAL069W 335 649
YAL068W-A 538 792
YAL067W-A 2480 2707
YAL066W 10092 10400
YAL064W-B 12047 12427
YAL064W 21526 21852
```

Answer: This requires rather a simple text manipulation trick. The key is that you need to know exactly what you want as a BED file. Read more about BED format [here](https://genome.ucsc.edu/FAQ/FAQformat.html#format1).
Question: Hello People,

I ran the new tuxedo pipeline for a reference based transcriptome analysis on HumanGenome19. After running HISAT2, I executed StringTie with -A option and got the abundance file and it looks like the following

Gene ID	Gene Name	Reference	Strand	Start	End	Coverage	FPKM	TPM

ENSG00000223972	DDX11L1	1	+	11869	14412	0.261954	0.056163	0.152539

ENSG00000227232	WASH7P	1	-	14363	29806	23.156389	5.585505	15.170202

ENSG00000243485	MIR1302-10	1	+	29554	31109	0.093882	0.048981	0.133031

Now what I would like to know is how is the "Coverage" value in 7th column calculated and in what units (or is it in percentage) ??

In the manual page, it shows "cov: The average per-base coverage for the transcript or exon" but still the calculation is not clear to me.

Answer: The coverage value shown in the output of StringTie (and in other genomics programs) is an average of all per-base coverages across the length of genomic segment (exon) or set of segments (transcript). For a gene (which is what is shown in the -A table), the union of all the exonic segments in the gene is considered for this calculation.

For example, an exon coverage value is calculated like this: all the read alignments intersecting the exon are considered and the coverage values for each base (i.e. the number of the read alignments covering that base) are added up and then divided by the length of that exon.

Additionally in StringTie there are a few other factors influencing that per-base coverage value (i.e. the way the number of alignments covering a base is actually "counted"):

- potential filtering of some of the read alignments (those considered suspicious/unreliable; they may appear in IGV but discarded by StringTie)
- weighing down the per-base coverage contribution of multi-mapped reads (i.e. if a read is mapped in *n* other places, we count its base coverage contribution as *1/n* instead of 1, for each "covered" base)
- for multi-transcript genes: distributing read alignments among overlapping transcripts according to the maximum flow algorithm

Note that the last factor above, although the most complex, does not influence *gene coverage* values, because at that level we do not need to worry about the distribution of read alignments between transcripts/exons (as mentioned above, simply considering the segment union of all exons in the gene is enough).
Question: Hi,

I'm trying to reproduce the calculation of moderated z-score (MODZ) described in PMID: 29195078 for some microarray data. Calculation is described as follows:

> Replicate-consensus signatures (MODZ) L1000 experiments are typically
> done in 3 or more biological replicates. We derive a consensus
> replicate signature by applying the moderated z-score (MODZ) procedure
> as follows. First, a pairwise Spearman correlation matrix is computed
> between the replicate signatures in the space of landmark genes with
> trivial self-correlations being ignored (set to 0). Then, weights for
> each replicate are computed as the sum of its correlations to the
> other replicates, normalized such that all weights sum to 1. Finally,
> the consensus signature is given by the linear combination of the
> replicate signatures with the coefficients set to the weights. This
> procedure serves to mitigate the effects of uncorrelated or outlier
> replicates, and can be thought of as a ‘de-noised’ representation of
> the given experiment’s transcriptional consequences.

I already calculated robust z-score for each gene (n=900) at every sample (n=6, 3 controls, 3 treatments) and I used the transposed matrix (genes as columns, samples as rows) to calculate pairwise spearman correlation using:

cor<- cor(tm, use="pairwise.complete.obs","spearman")

but I'm sure I'm missing something because I get error

Error in cor(tm, use = "pairwise.complete.obs", "spearman") :
'y' must be numeric

I don't know how to define x or y (I'm quite newby with R), but as I understand, I have to correlate control 1 to control 2, control 2 to control 3 and control 1 to control 3 individually and same for treatments and then compute the weights.

Any ideas about the code or how to proceed will be very much welcome!!!

Gema

Answer: As just mentioned on your other post, you need `method="spearman"`, since otherwise R is seeing the unnamed argument as the thing you should be correlating `tm` with. Also, you want `t(tm)`, as Kevin mentioned.
Question: Hi,

I have got the below files from the sequencing laboratory

```
Sample_L001_R1_001.fastq
Sample_L001_R1_002.fastq
Sample_L001_R1_003.fastq
Sample_L001_R1_004.fastq
Sample_L001_R1_005.fastq
Sample_L001_R1_006.fastq
Sample_L001_R1_007.fastq
Sample_L001_R2_001.fastq
Sample_L001_R2_002.fastq
Sample_L001_R2_003.fastq
Sample_L001_R2_004.fastq
Sample_L001_R2_005.fastq
Sample_L001_R2_006.fastq
Sample_L001_R2_007.fastq
Sample_L002_R1_001.fastq
Sample_L002_R1_002.fastq
Sample_L002_R1_003.fastq
Sample_L002_R1_004.fastq
Sample_L002_R1_005.fastq
Sample_L002_R1_006.fastq
Sample_L002_R1_007.fastq
Sample_L002_R2_001.fastq
Sample_L002_R2_002.fastq
Sample_L002_R2_003.fastq
Sample_L002_R2_004.fastq
Sample_L002_R2_005.fastq
Sample_L002_R2_006.fastq
Sample_L002_R2_007.fastq
```

I could understand that

`L**` represents lane, what else information about the sequencing can be drawn from these naming conventions. Could someone throw light?

Many Thanks

Answer: `Sample` is the sample name, `L???` is the lane, and `R?` is the read number (1 or 2), which is only meaningful for paired-end sequencing. The last 3 digits don't have any real meaning; the reads are just put into multiple smaller files instead of a single larger one. Note that, for example, `Sample1_L001_R1_001.fastq` and `Sample1_L001_R2_001.fastq` describe the same read pairs, so make sure to always process them together. The number of files per sample is then just determined by the sequencing depth and how large those files are set to be.
Question: I'm really liking RTG's vcfeval command, which takes in a gold standard vcf, a query vcf and a reference and outputs true positive vcf, true negative vcf, and the same for false positive, and false negative. My question is: is there some way to easily annotate the query vcf file so that each variant is annotated as TP, TN, FP and FN?

Answer: Yes, that's easy, just use `--output-mode annotate` (also check out the other output modes in case something else is more suitable to your needs).
Question: I aligned some paired-end RNA-seq data using STAR and then used featureCounts on the output .bam to assign them to genes. However I'm really confused as the output of featureCounts shows that I have more fragments (read-pairs) than the number of total input reads (I'm assuming also read-pairs/fragments) in the STAR output.

STAR output:

Number of input reads |	45661688
Average input read length |	199
UNIQUE READS:
Uniquely mapped reads number |	41983291
Uniquely mapped reads % |	91.94%

featureCounts output:

|| Paired-end reads are included. ||
|| Assign fragments (read pairs) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Read re-ordering is performed. ||
|| ||
|| Total fragments : 49756690 ||
|| Successfully assigned fragments : 36826889 (74.0%) ||

45661688 vs 49756690... where have these extra ~4 million reads come from? Also, aren't multimapped reads excluded from the .bam STAR output anyway so really it is 41983291 vs 49756690? Anyone any ideas?? Thanks!

Answer: You are mixing two things. STAR's "Log.Final.out" reports mapping statistics in reference to the number of input reads in the fastq file. So **Number of input reads | 45661688** refers to the total number of reads from the fastq file. Down the log file, STAR outputs several statistics which account for every input read:

- Uniquely mapped reads number | 41983291 (which you showed)
- Number of reads mapped to multiple loci
- Number of reads mapped to too many loci
- Number of reads unmapped: too many mismatches
- Number of reads unmapped: too short
- Number of reads unmapped: other
- Number of chimeric reads

Summing all these, you should recover the total number of input reads, that is, **the number of reads from the fastq file**.

featureCount log reports statistics in reference to the **number of mapped reads**, which can be different from the number of input reads from the fastq file for two reasons:

- unmapped reads may be left out of the bam file
- multi-mapping reads will increase the number of "fragments" reported by featureCounts

To get the total number of fragments, featureCounts counts multi-mapped reads. However, with default settings, featureCounts doesn't assign those counts to any feature - that is, multi-mappers results in zero counts to the features they map.

Regarding the questions from your comment:

> Would you recommend leaving the STAR default as 10 for aligning to the genome?

If you are using featureCounts for read summarization, this setting doesn't make any difference, as featureCounts will assign zero counts for multi-mappers. So, yeas, leave it at the default.

> On a side note I've also been reading a bit about mmquant, would you recommend using that instead of feature counts for assigning features?

No, I would stick to using featureCounts, or use a method which accounts for multi-mappers with a proper statistical model, such as RSEM, Salmon or kallisto.
Question: I am reading a paper that has a passage describing the pre-processing of gene expression data before conducting the experiments. The passage states "After conversion to a base-2 logarithm with a pseudo count of 0.125, batch normalization using ComBat was applied".

What exactly is a pseudo count? What I understood initially was that you add 0.125 to every value in your gene expression matrix and then take the logarithm of that to avoid taking the logarithm of 0 (which is not defined). This is based on my intuition though and I would like to know if this is correct and if there are other reasons why pseudo counts are used.

Answer: Your understanding is correct.

I personally like `log2(x+1)` because a pseudocount of 1 means you don't have to deal with negative numbers.
Question: I'm just biologist and not an expert in Bayesian statistics, but try to learn how to properly and wisely use the bayesian approach.
So, for my concatenated matrix (mtDNA, nuclearDNA + gaps in one of nuclear genes) I made partitioned analysis.

Fully resolved tree topology were yielded with matrix dissected on 9 partitions. To all of them were selected models of nucleotide substitutions (with jModelTest).

gene 1 1st nucleotide in codon

gene 1 2nd nucleotide in codon

gene 1 3st nucleotide in codon

gene 2 1st nucleotide in codon

gene 2 2nd nucleotide in codon

gene 2 3rd nucleotide in codon

gene 3 (rDNA)

gene 4 (rDNA)

gaps (as "standard" characters)

then

lset applyto=(1) nst=1 rates=gamma;
lset applyto=(2) nst=2 rates=equal;
lset applyto=(3) nst=6 rates=gamma;
lset applyto=(4) nst=1 rates=equal;
lset applyto=(5) nst=1 rates=propinv;
lset applyto=(6) nst=1 rates=propinv;
lset applyto=(7) nst=1 rates=gamma;
lset applyto=(8) nst=2 rates=equal;

After a while, I found Jeremy M. Brown's and Fredrik Ronquist's presentations, where they recommended not to choose models, but
"let the [bayesian] analysis sample different models... (reversible jump)"
and
"If you use ModelTest or MrModelTest: Do not fix parameters in MrBayes"

Does it mean for me, that I did everything wrong? With such too detailed partitition I reсieve distorted tree topology, I guess.
And can someone explain why is it need to find models by the Bayesian MCMC analysis itself?

Thank you for your attention. I will be glad to hear any comments regarding my analysis.

Answer: Hi Alice,

First of all, I do not think you did everything (or anything for that matter) wrong in the analysis.

Even with the "lset nst= mixed" command you still need to specify different models of rate variation across sites. So I think it is ok to use ModelTest models in a Bayesian analysis. With the obvious caveat that a model is just that, and can over- or underestimate what is really happening.

As to why you would want MrBayes to test all the models was suggested in their paper (Huelsenbeck et al., 2004):
"Current implementations of any of these criteria suffer from the limitation that only a small set of models are examined, or that the test does not allow easy comparison of non-nested models. In this article, we expand the pool of candidate substitution models to include all possible time-reversible models. This set includes seven models that have already been described."

I obviously do not know what tree you are running, but use common sense. The eyeball test (i.e. are things that should be together together ...) can go a long way.

Most importantly, make sure your trees have converged. I like to use AWTY (http://king2.scs.fsu.edu/CEBProjects/awty/awty_start.php) for that.
Question: Hi all, is it possible to determine the name of the reference file used from any VCF file?

Thanks in advance.

Answer: yes if there are some contigs lines in the VCF header.

##contig=

yes if there is a reference line in the header:

##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta

yes if you can get a rs# in the ID column with a column CHROM and POS, without any ambiguity

yes if you can check CHROM and POS vs a reference, without any ambiguity
Question: I've downloaded samfiles using conda
For converting sam to bam I used the command

samtools view -S -b /Users/dhee/Desktop/SRR12304924.s > /Users/dhee/Desktop/SRR12304924.bam

In getting an error as

dyld: Library not loaded: @rpath/libcrypto.1.0.0.dylib
Referenced from: /Users/dhee/miniconda3/bin/samtools
Reason: image not found
zsh: abort samtools view -S -b /Users/dhee/Desktop/SRR12304924.s

Answer: Try this and let me know if it fixes the issue:

conda uninstall samtools
conda update --all
conda install samtools
Question: I have several columns of gene symbols. I would like to know which gene exists similarly in each of columns

For example, I have three columns which shows that the first gene is similar in the first column and second column,

Is there anyone know how to check for this?

```
ADORA2B ADORA2B KLC1
AGPAT5 HOPX LEPR
ASS1 IGFBP7 LTBP3
C1QBP INHBA MBNL2
C4orf19 ITGB5 MLLT11
CASP1 ITGBL1 NOTCH3
CASP1 ITGBL1 NPR3
CASP1 NUAK1
CASP1 OLR1
CCL20 PDGFC
KBTBD11 PLA2G16
KLF4
ME2
MPDU1
NAT1
PBK
PSMB10
PSMB8
PSMB9
```

Answer: comm -1 -2 <(cut -f1 yourFile | sort) <(cut -f2 yourFile | sort)

Would compare the first and second columns and print only values that were in both columns.

comm -1 -2 <(cut -f1 yourFile | sort) <(cut -f3 yourFile | sort)

Would compare the first and third columns and print only values that were in both columns.

etc.

As far as I recall, by default sort uses dictionary order, which should work fine for your data. Or maybe the order they're in already matters? If yes, you should probably use awk. For example:

awk 'FS="\t" {if($1 == $2 && $2 == $3) print}' yourFile

Would only print lines where columns 1 2 and 3 have the same value
Question: One thing I really like about CWL is the ability to load CWL files into a CommandLineTool or Workflow object after generating the classes using `schema-salad-tool --codegen=python CommonWorkflowLanguage.yml > cwl_classes.py`. Is there a schema file that describes CWL inputs/job files similar to CommonWorkflowLanguage.yml? I have been working with these files as dicts, but would be super nice to be able to load them straight into a class object.

Answer: To expand a bit on what Michael said, the "inputs" and "outputs" section of every tool or workflow is a schema for the input object, so (although I have not tried it) it probably is not much more complicated than dumping the inputs section and using the code generator on it.