cannin/biostars_qa huggingface dataset:
Question: I have a `data.frame` with 48 columns. the columns are samples in triplicates, the rows are genes. The values are expression values.

For each group of triplicates I would like to calculate the averaged expression per gene, resulting in a new `data.frame` with 16 columns and for each row the average of the three triplicates for each group.

can this be done with the `tidyverse` tools?

thanks

a small example table of 6x9 is here

dput(head(normCounts[,1:9]))
structure(c(0, 4.89997034943019, 2.4499851747151, 0, 46.5497183195869,
14.6999110482906, 0.998187766715749, 1.9963755334315, 0, 0.998187766715749,
55.898514936082, 7.98550213372599, 0, 1.57112407949228, 0, 1.57112407949228,
53.4182187027374, 4.71337223847683, 0, 1.25548317693578, 0, 0,
52.7302934313026, 10.0438654154862, 0, 0, 0, 0, 66.3962127189125,
23.2386744516194, 2.18533123780511, 3.27799685670766, 0, 0, 65.5599371341532,
9.83399057012298, 0, 0, 0, 0, 74.1086143860152, 18.9580176336318,
0, 0, 0, 0, 66.8826789069951, 13.376535781399, 0, 0, 0, 0, 50.7776960416371,
13.0791035258762), .Dim = c(6L, 9L), .Dimnames = list(c("ENSMUSG00000103147",
"ENSMUSG00000102269", "ENSMUSG00000096126", "ENSMUSG00000102735",
"ENSMUSG00000098104", "ENSMUSG00000102175"), c("Sample_1", "Sample_2",
"Sample_3", "Sample_4", "Sample_5", "Sample_6", "Sample_7", "Sample_8",
"Sample_9")))

Answer: It will be a lot simpler to first make a group key table.

library("tidyverse")

groups <- tibble(
sample=colnames(normCounts),
group=rep(seq(1, ncol(normCounts)/3), each=3)
)

> groups
# A tibble: 9 x 2
sample group

1 Sample_1 1
2 Sample_2 1
3 Sample_3 1
4 Sample_4 2
5 Sample_5 2
6 Sample_6 2
7 Sample_7 3
8 Sample_8 3
9 Sample_9 3

Then you can pivot your count data to long format, join the groups, and group_by/summarize to get the means.

mean_exp <- normCounts %>%
as_tibble(rownames="gene") %>%
pivot_longer(starts_with("Sample"), names_to="sample", values_to="counts") %>%
left_join(groups, by="sample") %>%
group_by(gene, group) %>%
summarize(mean_count=mean(counts))

> head(mean_exp)
# A tibble: 6 x 3
# Groups: gene [2]
gene group mean_count

1 ENSMUSG00000096126 1 0.817
2 ENSMUSG00000096126 2 0
3 ENSMUSG00000096126 3 0
4 ENSMUSG00000098104 1 52.0
5 ENSMUSG00000098104 2 61.6
6 ENSMUSG00000098104 3 63.9

You can pivot back to a wider format if you want too.

mean_exp_wider <- pivot_wider(mean_exp, names_from=group, values_from=mean_count)

> mean_exp_wider
# A tibble: 6 x 4
# Groups: gene [6]
gene `1` `2` `3`

1 ENSMUSG00000096126 0.817 0 0
2 ENSMUSG00000098104 52.0 61.6 63.9
3 ENSMUSG00000102175 9.13 14.4 15.1
4 ENSMUSG00000102269 2.82 1.51 0
5 ENSMUSG00000102735 0.856 0 0
6 ENSMUSG00000103147 0.333 0.728 0
Question: Hi all,

I was wondering if the order of records in a bamfile produced by bwa aln/mem/sampe is guaranteed to be the same as in the fastq files that were used as input.

I checked the bwa manual but the only thing I found is this: "Repetitive read pairs will be placed randomly" (part of sampe description).

If anyone has an idea, I'd welcome your feedback!

Best,
~Lina

Answer: Yes, the same.
Question: I have a list of about 800 KEGG IDs for my genes of interest. I want to plot those genes and get a nice KEGG mapped plot for each of my 800 genes. How can I do this in R package? Thanks

Answer: You may use the `pathview` package:

1) use `kegg.gsets()` to download the KEGG pathways of your species.

2) map your gene ids to the pathways - I use `ids2indices()` from limma.

3) use `pathview()` to download KEGG pathways figures and color your genes on them.

Example:

library(pathview)
test.uid <- c("5494737" , "5495078", "5495093", "5494418", "5495039")
pv.out <- pathview(gene.data = test.uid, pathway.id = "ssl00190", species = "ssl",
out.suffix = "kegg.get.all", kegg.native = T)
Question: Hi,

I am a beginner of the RNA seq. Now I do some analysis about the paired-end data. For filtering the raw data, I use the Trimmomatic to classify the data to two parts: unpaired reads and paired reads.

So, I don't understand what's the unpaired. Because every reads should have its partner, I always get the same number of the Reads 1 and Reads 2. Why exist the unpaired reads?

Thanks for your reply.
Yunlong

Answer: When a trimming program trims data one of the reads (assuming you are using paired-end data) may become short and fail a criteria you have set (e.g. minimum length 25 bp). At the point that read is removed by the trimming program. Let us say that was from R2 file. A trimming program should remove the corresponding R1 read from R1-file (even though that may have passed) when R2 read is dropped. If it does not do that you are left with an unpaired read in R1 file.
Aligners expect reads to be in the same order in R1/R2 files. If they are not then you can get strange results (e.g. discordant mapping). Presence of unpaired reads in main sequence files (As @mastal points out trimmomatic should collect them in separate files) generally signifies improper use of a trimming program (or using a trimmer that is not paired-end aware). If that happens `repair.sh` tool from BBMap suite can be used to remove those unpaired reads and bring R1/R2 files back in sync.
Question: HI, im having trouble with this task
Write a perl script that will generate a new output file (“task1 output.txt”) which contains
the sequence name, length, and GC-content for each sequence. There should be a header
line which identifies the contents of columns (so the first line in the output file should be
“SeqName Length GC-Content” or something similar). The GC-content of a sequence is
defined as the percentage of bases that are G or C (from 0% to 100%), and a high GCcontent
is associated with coding sequences.

>Seq1
ACGT
Then your output file should look like:
SeqName Length GC-Content
Seq1 4 50

I can do the in and out for the file handles, but im confused as to what to put in my while loop. and how will it know to match in the file ?

Answer: Like mentioned, GC content is the percentage of bases that are G or C in the sequences. Percentage is calculated quite easily using basic math once you obtain the counts of bases that are G or C in and the total sequence length for each sequence. Your loop with iterate through the file, executing its code block for each sequence it finds.

I cannot - and I hope others do not too - provide code. That would cripple learning.
Question: Hi everybody,

I have a list of recurrent mutations from a WGS dataset. I would like to know on which exons they sit (the exon number for a certain transcript).

If there are multiple transcripts at a certain position, annovar puts out a list with all nonsynonymous SNVs and I can just pick the transcript I want. My only problem is the following: If there is for example a stoploss for any of the transcripts at that genomic position, ANNOVAR doesn't output the nonsynonymous SNVs anymore, but only a list of mutation of that category, which has a higher precedence.

Is there any way to change that? Can I either change the precedence, or tell the software to output the effects on all transcripts or anything similar?

Best regards, and thanks a lot in advance,
Gero

Answer: Hello, all that you need to add is the `--separate` command line parameter:

#1, view input:

cat test.ann
3	38182727	38182727	A	C
3	38182316	38182316	A	G
3	38182641	38182641	T	C

#2, annotate:

perl annotate_variation.pl -out ex1 -build hg19 test.ann /Programs/annovar/humandb/ --separate

#3, view output:

cat ex1.exonic_variant_function

line1	nonsynonymous SNV	MYD88:NM_001172568:exon4:c.A745C:p.T249P,MYD88:NM_002468:exon5:c.A880C:p.T294P,MYD88:NM_001172567:exon5:c.A904C:p.T302P,	3	38182727	38182727	A	C
line2	nonsynonymous SNV	MYD88:NM_001172568:exon3:c.A617G:p.K206R,MYD88:NM_001172569:exon3:c.A571G:p.N191D,MYD88:NM_002468:exon4:c.A752G:p.K251R,MYD88:NM_001172566:exon2:c.A436G:p.N146D,MYD88:NM_001172567:exon4:c.A776G:p.K259R,	3	38182316	38182316	A	G
line3	stoploss	MYD88:NM_001172569:exon4:c.T613C:p.X205R,MYD88:NM_001172566:exon3:c.T478C:p.X160R,	3	38182641	38182641	TC
line3	nonsynonymous SNV	MYD88:NM_001172568:exon4:c.T659C:p.L220P,MYD88:NM_002468:exon5:c.T794C:p.L265P,MYD88:NM_001172567:exon5:c.T818C:p.L273P,	3	38182641	38182641	T	C
Question: Hi!

I'm starting to use kallisto to do transcript-level expression quantification.
I have some questions:

**1)** Does `kallisto` infer the strandness of the input data just like `salmon` does (`--libType A`)? I guess the answer is no.

**2)** For other hand, `kallisto` has the next to options:

--fr-stranded Strand specific reads, first read forward
--rf-stranded Strand specific reads, first read reverse

Are these options only working for PE data?

**3)** Regarding the fragment length estimation when using SE datasets:

-l, --fragment-length=DOUBLE Estimated average fragment length
-s, --sd=DOUBLE Estimated standard deviation of fragment length
(default: -l, -s values are estimated from paired
end data, but are required when using --single)

What does `DOUBLE` mean? Do we have to specify the double of the number calculated?

Thank you in advance

Answer: 1. No
2. They should work for SE data too (never tried, though). You probably want --rf-stranded for anything remotely recent.
3. An example of a double is `200.0` or `123.4`. That is, any number with a decimal point. The documentation there should really be changed, since I don't expect those without C/C++/etc. programming experience to know that "double" means "double precision floating point value" (or what that even means)).
Question: Who can give a finally determined answer: **what is paired and unpaired reads by trimmomatic**? what kind of sequences are in unpaired reads? please see [the question link and my comments][1]. Many many many thanks! here is my comments:

> "Hi Genomax, I see your reply for this question, but I still do not
> understand what is unpaired reads or unpaired.fastq file? based on
> your answer, my understanding is that for a Paired End sequencing,
> generally, the types of sequences in the R1 file is equal to that in
> the R2 file. Here, we do not care about the number of each sequence.
> for instance, if one sequence cannot pass the QC (set in
> trimmomatic)in R1 file, but this sequence pass the QC in R2 file,
> however, this sequence in both R1 and R2 file will be classified into
> unpaired reads/.fastq file, which means all the copies in R1 and R2
> files also will be classified into the unpaired reads. Or another
> understanding is that one sequence exist in both R1 and R2 file, but
> one copy in either R1 or R2 cannot pass the QC, this copy will be
> classified into unpaired fastq file/reads. (I think the second view
> might right). if so, some guys also mentioned using the unpaired for
> alignment/mapping with BWA, whether these under-QC sequences should be
> dealt with trimmomatic again with a strict set? are they useful?
> Finally, what kind of sequences are in unpaired.fastq file? please
> five some examples. my email is wangqk198738@163.com Thanks a lot!"

[1]: https://www.biostars.org/p/192788/

Answer: > my understanding is that for a Paired End sequencing, generally, the
> types of sequences in the R1 file is equal to that in the R2 file.

I don't know what you mean by that. There should only be one type/format of sequences, fastq. Number of reads in `R1/R2` files will be identical when the sequence comes off the sequencer.

> Here, we do not care about the number of each sequence.

Some may not (if you are willing to accept files that are out of sync in terms of order of `R1/R2` reads. Aligners will produce odd results if you use such a file.

> if one sequence cannot pass the QC (set in trimmomatic)in R1 file, but
> this sequence pass the QC in R2 file, however, this sequence in both
> R1 and R2 file will be classified into unpaired reads/.fastq file,
> which means all the copies in R1 and R2 files also will be classified
> into the unpaired reads.

That is the reason you should trim reads together and use a trimming program that is PE aware. There should be only one copy of sequence for each cluster coordinate.
Question: I recently used BLAT to generate a large number of alignments but the input files were on the chromosome level, not individual genes. I realize BLAT and filtering the subject by taxonomical ID is an option, but BLAT seems to be a stricter option. The alignments are 60bp with no gaps. What would be the best way to find the gene associated with each of these regions, if it exists?

Answer: As @AndreiR as said...

Download a bed file of all the genes of interest (from the ucsc table browser you can download ready-made bed files)

then make a bed file from your alignments:

then install bedtools and:

intersectBed -a alignments.bed -b refseq.bed -wo > intersections.txt

This should give you all the info you need.
Question: I have a biopython program that uses newick NHX extensions to pass information to the "R" tree plotting program ggtree. I am certain this used to work, but with the current version of biopython, it does not seem possible to write arbitrary comments in newick output.

(1) Am I confused, has BioPython phylo.write('newick') never supported NHX comments?
(2) is there an option for phylo.write() to write out the comment string?

Answer: This is occurring because of a typo in the BioPython Phylo/NewickIO.py code (line 87), where the word "comment" is mis-spelled "coment". Fixing the typo fixes the problem. I have opened an issue on GitHub.
Question: I just want to make it clear.

I need to calculate FPKM. I use this formula:
**Normalized = [(raw_read_count)(10^9)] / [(gene_length)(XXXX)],**

*XXXX = the count of all reads that are aligned to protein-coding genes in that alignment.*

How should I calculate XXXX? Is it just sum of **all** raw_read_counts after htseq-count (e.g. in R it will be XXXX <- sum(collumn_with_raw_red_counts)?

Thanks!

Answer: I've got answer from GDC portal:

> 1. Download GTF files used in HTSeq analyses: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files (GDC.h38 GENCODE v22 GTF)
2. Extract only protein-coding gene IDs: less gencode.v22.annotation.gtf | grep "\tgene\t" | grep protein_coding | cut -f9 | cut -f2 -d '"' > EnsembleIDsPCG.txt
3. Use resulting list to extract only protein-coding values from counts file: less CountFile.txt | grep -Ff ProteinCodingGeneList.txt > CountOnlyProt.txt
4. Sum the values of "CountOnlyProt.txt" and that will give you your denominator value.

My problem was that I counted reads for all genes, but should only for protein-coding.

P.S. thanks to GDC support team!
Question: Hi,

I need the .dat file for the UniProt Metazoa. .dat has the EMBL format, which I need. But I don't see it here:
**https://www.uniprot.org/taxonomy/33208#**

I can download the clade specific .dat files from here: **https://www.uniprot.org/downloads**. But I don't see the metazoan .dat file here.

Cab anybody help?

Cheers!!

Answer: Go to this page: https://www.uniprot.org/uniprot/?query=reviewed:yes%20taxonomy:33208 and click on the Download button. If you choose `Text` from the drop-down list for Format, you will get the data in EMBL format.
Question: Hello everyone

MuTect2 from GATK3 discarded `87.33%` of the reads. This is during preparing panel of normals.
I will use sp1.vcf in the PON creation !

sp1.vcf: sp1.bam
java -jar ${GATK}/GenomeAnalysisTK.jar \
-T MuTect2 \
-R ${hg38}.fasta \
-I:tumor sp1.bam \
--dbsnp ${DBSNP} \
--cosmic ${COSMIC} \
--artifact_detection_mode \
-o sp1.vcf

**Does this seem okay to you? Is there anything I can do to fix non primary reads `(42.77%)` and duplicate reads `(44.56%)`?**

INFO 00:57:17,188 MicroScheduler - 158885515 reads were filtered out during the traversal out of approximately 181931389 total reads (87.33%)
INFO 00:57:17,190 MicroScheduler - -> 9210 reads (0.01% of total) failing BadCigarFilter
INFO 00:57:17,192 MicroScheduler - -> 81066454 reads (44.56% of total) failing DuplicateReadFilter
INFO 00:57:17,193 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
INFO 00:57:17,195 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
INFO 00:57:17,197 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
INFO 00:57:17,198 MicroScheduler - -> 77809851 reads (42.77% of total) failing NotPrimaryAlignmentFilter
INFO 00:57:17,200 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter
------------------------------------------------------------------------------------------
Done. ------------------------------------------------------------------------------------------

**This is how STAR log summary looks:**

UNIQUE READS:
Uniquely mapped reads number | 29463720
Uniquely mapped reads % | 72.35%
Average mapped length | 271.00
Number of splices: Total | 21158459
Number of splices: Annotated (sjdb) | 21067563
Number of splices: GT/AG | 20952149
Number of splices: GC/AG | 126813
Number of splices: AT/AC | 5622
Number of splices: Non-canonical | 73875
Mismatch rate per base, % | 0.48%
Deletion rate per base | 0.02%
Deletion average length | 1.27
Insertion rate per base | 0.01%
Insertion average length | 1.84
MULTI-MAPPING READS:
Number of reads mapped to multiple loci | 9507878
% of reads mapped to multiple loci | 23.35%
Number of reads mapped to too many loci | 355263
% of reads mapped to too many loci | 0.87%
UNMAPPED READS:
% of reads unmapped: too many mismatches | 0.00%
% of reads unmapped: too short | 2.76%
% of reads unmapped: other | 0.67%
CHIMERIC READS:
Number of chimeric reads | 0
% of chimeric reads | 0.00%

Answer: If you are sure you want to retain non primary and duplicates, I think you can add to mutect the option `--disable_read_filter NotPrimaryAlignmentFilter --disable_read_filter DuplicateReadFilter` (not tested). See some docs [here.][1]

[1]: https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_cancer_m2_MuTect2.php
Question: Hi,

We have experimentally found that given gene X causes some phenotype. To start exploring the underlying mechanism, we decided to use TCGA data as its freely available. So within the cancer patients can we simple divide two samples based on gene X expression.

The first thing that came to my mind was Z normalise X gene and get the [-inf, +1.96] and [196,+inf] patients and do differentially expressed gene analysis.

Pleas keep in mind that this is only for preliminary data. So we aim just to hand wave and explore.
Any idea would make my day? Please.

Best regards,

Tunc.

Answer: Do I understand correctly that you are dividing the samples by the expression of the gene (rather than observing the expression of the gene in different groups)?

If so, the first thing I would do is make a histogram of the log TPM and see if it is bimodal. Then if it is you will know where the cutoff is. If it's not bimodal it may be just a noisy gene. I would do it in GTEX for the the normal tissue as well.
Question: I used The 9th column of `.bam` file to get the insert size of atac-seq data.

```
samtools view bam/${FILE}"_mm10.bam" | awk '$9>0' | cut -f 9 | sort | uniq -c | sort -b -k2,2n > insert_size_plotting/${FILE}".txt"
```

But I found all of the inserts are less than 500bp,
From the insert size plot, I can see there should be longer than 500bp inserts, but why bowtie2 doesn’t report >500bp inserts?
![enter image description here][1]
Thanks

[1]: /media/images/c6891381-0797-4efd-8f9f-5d375f85

Answer: By default, bowtie2 expects a maximum fragment length of 500 bp for paired-end alignments. To adjust this, change the `-X` or `--maxins` parameter.

http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: R Programming for Bioinformatics (Chapman & Hall/CRC Computer Science & Data Analysis)

Bioinformatics and Computational Biology Solutions Using R and Bioconductor (Statistics for Biology and Health)

Quick-R

cran: TheRGuide

cran: An Introduction to R

Google's R Style Guide

A Beginner's Guide to R

R Tutorial Series: R Beginner’s Guide and R Bloggers Updates
Question: I have a list of gene name in a file

```
CHRNB2
EGR2
GCK
KRT14
LMNA
FGF3
TK2
ABCC8
```

How can U map them to uniprot ID?

**P.S** I tried Uniprot "ID mapping" (from-"GENEID" to-"UNIPROTKB AC").But it couldn't map.

Please suggest me what to do.Thnx

Answer: Use Mygene.info. You can do batch requests via post, or you can use the live API to do batch requests as well.

Here's how via the live API:

Click on the "Try API live!", select "gene query service". Click on "post"

For "q" put in your gene names separated by a comma.

For "scopes" type "symbol" (without the quotation marks)

For "fields", use "symbol,entrezgene,uniprot" and any other parameter of interest

Click "try it" when done.

Result will be in the response body.
Question: Can anyone share me a good GATK4 workflow, I Couldnt find anything good that explanins the process in a crisp and easy way.

Answer: Completly, agree - there's a real lack of clarity when it comes to documentation of some tools and options. And it's only trolling through dozens of forum posts on any number of issues that you can begin to troubleshoot common issues. Example workflows won't resolve those headaches but the following resources may help you get started.

1. IBM White Paper - I found this most useful because I didn't want a .wdl cromwell workflow. https://www.ibm.com/downloads/cas/ZJQD0QAL
2. NIH HPC has a Good GATK4 tutorial/workflow with detailed explanations and examples I wish had been available when I made my workflow - https://hpc.nih.gov/training/gatk_tutorial/
3. Broad github .wdl scripts - Not ideal if you don't want to use .wdl, but they can still clarify some things. Making them even less ideal, the Broad regularly deems these scripts as 'outdated'. https://github.com/gatk-workflows/broad-prod-wgs-germline-snps-indels and also https://github.com/gatk-workflows/gatk4-basic-joint-genotyping/blob/main/gatk4-basic-joint-genotyping.wdl
4. Workflow using Nextflow for the variant calling - https://github.com/gencorefacility/variant-calling-pipeline-gatk4