cannin/biostars_qa huggingface dataset:
Question: Hi there,

I've used VEP (GRCh38.90) with the option --canonical to annotate my variants from WES and want to focus on canonical transcripts. This gives me a column if the associated transcript is the canonical one.

Now I'd like to retrieve a list of all canonical transcripts as used by VEP and failed to do so. I've tried to get all canonical transcripts following this thread https://groups.google.com/forum/#!topic/biomart-users/skO4zgqzGBA and I've downloaded knownCanonical.txt from UCSC.

However, both approaches have discrepancies to VEP. For example VEP reports ENST00000288602 as the canonical transcript of BRAF, while the output from biomart reports ENST00000496384 and knownCanonical reports ENST00000646891. I would like to go with VEP, so is there a way to obtain such a list from the VEP cache for example?

Any help is greatly appreciated.

Answer: [Glossary](http://www.ensembl.org/info/website/glossary.html)

**Canonical transcript**
The canonical transcript is used in the gene tree analysis in Ensembl and does not necessarily reflect the most biologically relevant transcript of a gene. For human, the canonical transcript for a gene is set according to the following hierarchy:

1. Longest CCDS translation with no stop codons.

2. If no (1), choose the longest Ensembl/Havana merged translation with no stop codons.

3. If no (2), choose the longest translation with no stop codons.

4. If no translation, choose the longest non-protein-coding transcript.
Question: Consider the following example from [bedtools intersect help page][1]. Say, there are 4 bed files as follows

$ cat query.bed
chr1 1 20
chr1 40 45
chr1 70 90
chr1 105 120
chr2 1 20
chr2 40 45
chr2 70 90
chr2 105 120
chr3 1 20
chr3 40 45
chr3 70 90
chr3 105 120
chr3 150 200
chr4 10 20

--

$ cat d1.bed
chr1 5 25
chr1 65 75
chr1 95 100
chr2 5 25
chr2 65 75
chr2 95 100
chr3 5 25
chr3 65 75
chr3 95 100

--

$ cat d2.bed
chr1 40 50
chr1 110 125
chr2 40 50
chr2 110 125
chr3 40 50
chr3 110 125

--

$ cat d3.bed
chr1 85 115
chr2 85 115
chr3 85 115

If we use bedtools like this

$ bedtools intersect -a query.bed \
-b d1.bed d2.bed d3.bed
chr1 5 20
chr1 40 45
chr1 70 75
chr1 85 90
chr1 110 120
chr1 105 115
chr2 5 20
chr2 40 45
chr2 70 75
chr2 85 90
chr2 110 120
chr2 105 115
chr3 5 20
chr3 40 45
chr3 70 75
chr3 85 90
chr3 110 120
chr3 105 115

It reports the interval for e.g. `chr1 5 20` , however, I don't want that to be reported, because it is not present in `d3.bed`. I know that it is being reported because it is overlapping in `query.bed` and `d1.bed`. But my requirement is different. How can I tweak `bedtools` to get ONLY those records which are overlapping across all files.

[1]: https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

Answer: Try the following

bedtools multiinter -i *.bed > test1
files="$(ls *.bed|wc -l)"
awk -v OFS='\t' -v files=$files '$4==files' test1 |cut -f 1-3 > test2

What it does:

bedtools multiinter -i *.bed > test1 # use bedtools multi intersect to get statistics on intersections of all bed files
files="$(ls *.bed|wc -l)" # make a variable called files that is the number of bed files in the directory
awk -v OFS='\t' -v files=$files '$4==files' test1 |cut -f 1-3 > test2 # get only regions that are found in all of the files (column 4 is num, which is the number of input files that have that region)

test2 is the file that contains regions found in all files assuming that the only bed files in the directory are query.bed, d1.bed, d2.bed, and d3.bed
Question: Hello

I have a few of questions while studying tool MARKDUPLICATES metrics file.

A list of `SECONDARY_OR_SUPPLEMENTARY_RDS` is there.

I know SECONDARY READS, but I am wondering about SUPPLEMENTARY READS.

Thanks in advance for full explanation.

Answer: Section 1.2 "Terminologies and Concepts" of the [SAM file format specifications][1] explain this:

> Chimeric alignment
>
>An alignment of a read that cannot be represented
> as a linear alignment. A chimeric alignment is represented as a set of
> linear alignments that do not have large overlaps. Typically, one of
> the linear alignments in a chimeric alignment is considered the
> “representative” alignment, and the others are called “supplementary”
> and are distinguished by the supplementary alignment flag. All the SAM
> records in a chimeric alignment have the same QNAME and the same
> values for 0x40 and 0x80 flags (see Section 1.4). The decision
> regarding which linear alignment is representative is arbitrary.

[1]: https://samtools.github.io/hts-specs/SAMv1.pdf
Question: Hi all,

I am using clusterProfiler to perform KEGG enrichment. I have a list of gene Symbol (Rat genes). First I need to translate the symbols to EntrezID.

I installed the Genome wide annotation for Rat:

source("https://bioconductor.org/biocLite.R")
biocLite("org.Rn.eg.db")

Then I tried to run the command:

x <- c("GPX3", "GLRX", "LBP", "CRYAB", "DEFB1", "HCLS1", "SOD2", "HSPA2",
"ORM1", "IGFBP1", "PTHLH", "GPC3", "IGFBP3","TOB1", "MITF", "NDRG1",
"NR1H4", "FGFR3", "PVR", "IL6", "PTPRM", "ERBB2", "NID2", "LAMB1",
"COMP", "PLS3", "MCAM", "SPP1", "LAMC1", "COL4A2", "COL4A1", "MYOC",
"ANXA4", "TFPI2", "CST6", "SLPI", "TIMP2", "CPM", "GGT1", "NNMT",
"MAL", "EEF1A2", "HGD", "TCN2", "CDA", "PCCA", "CRYM", "PDXK",
"STC1", "WARS", "HMOX1", "FXYD2", "RBP4", "SLC6A12", "KDELR3", "ITM2B")

eg = bitr(x, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Rn.eg.db")

Then I got the warning message:

> 'select()' returned 1:1 mapping between keys and columns Warning
> message: In bitr(x4, fromType = "SYMBOL", toType = "ENTREZID", OrgDb =
> Rat) :
> 98.21% of input gene IDs are fail to map...

Could anyone tell me what happened and how should I do?

Many thanks,
Stanley

Answer: The root of your problem is that your gene names are all upper case, but for rat (and mouse etc) generally only the first letter is capitalised, and the `bitr` conversion is case sensitive. We can use the function `str_to_title` in the **stringr** package to fix this:

install.packages('stringr')
x2 <- stringr::str_to_title(x)
head(x2)

Here's what they look like now:

[1] "Gpx3" "Glrx" "Lbp" "Cryab" "Defb1" "Hcls1"

Then rerun your query:

eg = bitr(x2, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Rn.eg.db")

and we get a much better proportion of converted IDs

'select()' returned 1:1 mapping between keys and columns
Warning message:
In bitr(str_to_title(x), fromType = "SYMBOL", toType = "ENTREZID", :
1.79% of input gene IDs are fail to map...

Looking at some of the other answers, **PVR** is an exception to the capitalisation rule. Since Kevin's **biomaRt** answer isn't case sensitive, it's perhaps the best solution assuming you have reliable internet.
Question: I am a complete newbie to machine learning and this is my first try, but I am not sure I am approaching this field in the right way so would very appreciate feedback.

My data and what I am trying to answer:
I have a data frame filled with IDs measured across three different features, a classification based on IDs, measured Features I would like to test, and values that depend on the particular ID-Feature combination.

My goal is to test whether one Feature, or a combination of Features, can be predictive of the Class_Name (defined from the ID). For now the data is in long format as below (example, the number of IDs varies for each feature and it is in the order of 100, so bigger than in this example):

df <- data.frame (Name = c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6", "ID9", "ID5", "ID7", "ID8", "ID9", "ID10", "ID4", "ID11", "ID12", "ID8", "ID9", "ID11", "ID13", "ID8", "ID12", "ID4"), Class_Name = c("orange", "orange", "red", "orange", "blue", "blue", "red", "blue", "orange", "blue", "red", "red", "orange", "orange", "red", "blue", "red", "orange", "red", "blue", "red", "orange"), Features = c("F1","F1","F1", "F1","F1","F1","F1","F2","F2","F2","F2","F2","F2","F2","F2","F3", "F3", "F3", "F3", "F3", "F3", "F3"), Values = c(21, 32, -36, 11, -62, 32, 34, -21, -43, -68, -24, -19, 28, 33, -33, 15, 2, 13, -99, -86, 3, 0))
df$Class_Name= as.factor(df$Class_Name)

I think I need to test each Feature in turn (script below). However, how do I then test whether combinations of features help with the prediction?

I am following tutorials on caret particularly [this][1]. (If anyone knows of tutorials in R with examples like this one, please post them), but I am afraid I am not approaching my question in the correct way, and I was wondering if you could direct me to the correct path.

Script I have tried:

library(caret)
# library(dplyr)
# table(df$Class_Name)
# select 2 so have balanced set for each class?
# df <- df %>% group_by(Features, Class_Name) %>% slice_sample(n=2) %>% data.frame() # this will be 100 in my data

control <- trainControl(method="cv", number=10)
metric <- "Accuracy"

features = unique(df$Features)
res = data.frame()
for (f in features) {
dataset = df[df$Features == f, ]

# Random Forest
fit.rf <- train(Class_Name~Values, data=dataset, method="rf", metric=metric, trControl=control)
print(fit.rf)
}

[1]: https://machinelearningmastery.com/machine-learning-in-r-step-by-step/

Answer: For testing multiple interactions, you can use:

model = randomForest(response ~ (var1 + var2 + var3)^2, data = your_data)
See a similar question from Stack Overflow here:
[Interaction Effects in R][1]

[1]: https://stackoverflow.com/questions/39886562/how-to-check-interaction-effects-for-a-lot-of-predictors-in-r
However, I believe caret is automatically looking for interaction effects.
To see which variables are contributing the most to your model, use varImp().
You can also check out the book "An Introduction to Statistical Learning" by James et al., which should answer all of your questions.
Question: I am trying to make UpSetR plot for nine vcf files using [working code][1]. However, the output plot shows only five files out of total nine files.

[1]: https://www.biostars.org/p/379679/#379721

Answer: Try increasing `nsets` when calling `upset`, like:

```
upset(fromList(read_sets), order.by = "freq", nsets = 9)
```

`?upset` shows:
>nsets	Number of sets to look at
Question: Dear all,

Do you have any idea how to easy extract contigs from fasta file wich contains specific sequence?

For example:

My sequence:

ACCGTACCC

My FASTA:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
>c1044
ATGAATAAAATAATTTTGTATCATAAATCGAGCTGTTAATTATT
>c1044
TTCATATTTGTAGCTAAGCAGAGGCGAAGCGTTCTTGTATCG
```

My output:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
```

Thank you so much for any ideas and help.

Answer: Use biopython or bioperl. With biopython, you could either use the re module or even just `find()` on the `str()` representation of each sequence. Either of these should be simple enough if you're familiar with either perl or python.
Question: I'm at ftp://ftp.ncbi.nih.gov/entrez/misc/data/gc.prt with all the genetic codes. Trying to figure out what genetic code *Nostoc punctiforme* uses. Is there a database that has which version they use? I know it's not standard and they have multiple start codons.

Answer: If you go to NCBI Taxonomy and search for Nostoc punctiforme, you'll reach its Taxonomy page with links for several strains, then clicking on the top level Nostoc punctiforme link will take you to the page which will have the genetic code you are looking for.
Question: I think I missed a detail in my understanding of missing residues. I understand that in some cases residues don't end up as part of the PDB model because of crystallization issues. However, aside from that case, will numbers be missing in the sequence of residues even though the residues are contiguous for the same reason that insert records may be used in the opposite scenario? In those cases, during processing they should be treated as contiguous?

Is this part of the reason why software packages such as Bio3d offer features for detecting chain breaks? Is this sort of software applicable to distinguishing the above cases?

Answer: Yes, as you wrote, the residue sequence number in the PDB format doesn't need to be sequential.
It's well explained by F.C. Bernstein here:
https://lists.sdsc.edu/pipermail/pdb-l/2004-March/001513.html
And some atoms or residues may be missing from the model.

How to distinguish the two cases?
If the PDB file has the SEQRES record - it should contain the list of all residues in the chain -- also the ones missing from the model. So this is one way of checking for missing residues.
Alternatively, you may check distances between atoms.
(I don't know anything about Bio3d).
Question: Hi guys,

I have protein expression data in a data frame df, where the proteins are rows and columns are sample ids with abundance values. Such as:

Sample 1 Sample 2 Sample 3
RPH3A
CA11
AIFM1

I want to keep those proteins that I have data for in at least 50 % of the samples. Any help would be appreciated?

Answer: You can calculate the percentage of NAs in each row and add that as a column to your data frame in R. Then you can use awk to remove rows that have more than 50% NAs.

In R:

count_na_func = function(x) {
sum(is.na(x))
}

df$na_percent = (apply(df, 1, count_na_func))/(ncol(df) - 1)
write.table(df, file = "datana.tsv", row.names=FALSE, sep="\t")

You can use some practice data to make sure that the na_percent column is computed correctly in R.

In bash:

awk '{if($12<=0.5){print}}' datana.tsv > newdata.tsv

Note that R is 1-indexed, while bash is 0-indexed. na_percent was in column 12 of my data frame, but it will probably be in a different column in yours, so you can substitute that column number for 12 in the awk code.
Question: Hi

I have a list of chromosomes and positions that looks like this:

1	10045
1	93056
1	109272
1	127711
1	127822
.
.
.

And now I would like to use it to remove them from my vcf file. Do you know how to do this?

Answer: [bcftools][1] can do this:

$ bcftools view -T ^list_snp_exclude.txt input.vcf > output.vf

With the `^` before the file with the coordinates one tell `bcftools` to exclude these regions.

fin swimmer

[1]: https://www.htslib.org/doc/bcftools.html#common_options
Question: Hi Everyone,

I am working on Whole Genome Sequencing and analysis of Human genome from illumina HiSeq platform with about 30X coverage. Each sample (human genome) have about 250-300 fastq.gz files, whom I am dealing with 'fastqc' for quality check using following command :

/usr/local/bin/fastqc -t 8 -f fastq -o OUT/ -casava *.gz -noextract

Although it is running fine and generating equal number of "fastqc.zip" files which I unzipped using unzip '*.zip'. So, here I have 2 questions:

1. Can I merge two or more fastq files and then run fastqc on those merged files? If yes, how should I merge those fastq files?
2. I have to manually check 250-300 fastqc folder to know the quality by opening .html page. Is there any way where I can have summary of overall quality of the fastq files in a flowcell?

Please let me know your comments. I'll be highly thankful to you.

Best,
Ravi

Answer: Not only can you merge the fastq files but your life might be easier if you do. For merging them, a simple `cat` will suffice. I should note that you don't have to be delivered 300 some odd files, you can request that whomever is doing the sequencing just give you a two files (assuming paired-end) per sample/library (the `bcl2fastq` program that they use to process the bcl files produced by the sequencer can trivially do this).

If you don't want to wait until all of the files are merged, you can likely just use a named pipe as input to fastqc. Something like:

mkfifo foo.fastq.gz
cat sample_L1_R1_???.fastq.gz > foo.fastq.gz
fastqc foo.fastq.gz

Given that fastqc is written in java, I can't guarantee that it'll properly handle block gzipped files like that (the java gzip library has been broken for years). You can always `zcat` instead. I should note that the only reason process substitution likely wouldn't work is that fastqc names the output files after the input file name.

For 2. it depends a bit on what you want. The sequencing facility actually has an idea about that already (it's produced by the machine). It's easy enough to just ask them (they can also give you a break down of how many reads per sample, their average quality (also per sample), etc.). For our internal pipeline, I have a pdf produced with that sort of information, since it's a bit quicker to look first at a single table like that than to trudge through all of the fastqc files. BTW, fastqc also produces an HTML file with the images included. When I QC flowcells before sending results to our local groups those are what I personally look at...it's quicker than dealing with the zip files.
Question: Hello,

I have a matrix in which I am trying to retrieve specific rows from it and save it in text

An example matrix is

```
Affy ID DDM1 DGKI2 FDGYYY1 GUHIL6
1438_at 0.0635 0.2065 -0.2112 0.0856
1487_at 0.071 -0.1315 0.0263 0.0198
1494_f_at 0.0045 -0.0237 0.0156 -0.1352
1598_g_at -0.0541 0.0006 -0.1369 -0.0589
160020_at -0.0925 0.2182 -0.1967 -0.0074
1729_at -0.0017 -0.2209 -0.086 -0.0709
1773_at -0.0273 -0.0181 0.1042 0.0136
177_at -0.0276 -0.2563 0.3975 -0.0535
179_at -0.0472 0.0979 -0.216 -0.2814
1861_at 0.0121 -0.4038 0.0016 0.0334
200000_s_at -0.1021 -0.0887 -0.0452 0.0035
```

Lets say the name of this matrix is `M.txt` and the selected rows is in a list named `mSelected.txt` (consisting of `1438_at`, `1729_at` and `200000_s_at`). My output should look like the following file

```
Affy ID DDM1 DGKI2 FDGYYY1 GUHIL6
1438_at 0.0635 0.2065 -0.2112 0.0856
1729_at -0.0017 -0.2209 -0.086 -0.0709
200000_s_at -0.1021 -0.0887 -0.0452 0.0035
```

Is there also anyway to convert their Affy ID to Gene name?

```
$ head /Users/Desktop/mSelected.txt | cat -vet
1438_at^M1729_at^M200000_s_at

$ head /Users/Desktop/m.txt | cat -vet
Affy ID ^I DDM1 ^IDGKI2 ^IFDGYYY1 ^I GUHIL6^M1438_at^I0.0635^I0.2065^I-0.2112^I0.0856^M1487_at^I0.071^I-0.1315^I0.0263^I0.0198^M1494_f_at^I0.0045^I-0.0237^I0.0156^I-0.1352^M1598_g_at^I-0.0541^I0.0006^I-0.1369^I-0.0589^M160020_at^I-0.0925^I0.2182^I-0.1967^I-0.0074^M1729_at^I-0.0017^I-0.2209^I-0.086^I-0.0709^M1773_at^I-0.0273^I-0.0181^I0.1042^I0.0136^M177_at^I-0.0276^I-0.2563^I0.3975^I-0.0535^M179_at^I-0.0472^I0.0979^I-0.216^I-0.2814^M1861_at^I0.0121^I-0.4038^I0.0016^I0.0334^M200000_s_at^I-0.1021^I-0.0887^I-0.0452^I0.0035
```

Answer: awk 'FNR==NR{a[$0];next}($1 in a)' mSelected.txt M.txt

If the files are huge, join could be the fastest way:

join -t "*" -1 1 -2 1 -o 1.1,1.2,1.3,1.4,1.5 <(sort -k1,1 M.txt) <(sort -k1,1 mSelected.txt) > output

*literal tab = ctrl-v-tab
Question: I am running bowtie2 for mate pairs specifying this parameter:

bowtie2 --local --phred64 --threads 25 -x "rest of the command"

It gives me this error:

```
Saw ASCII character -126 but expected 64-based Phred qual.
Try not specifying --solexa1.3-quals/--phred64-quals.
terminate called after throwing an instance of 'int'
bowtie2-align died with signal 6 (ABRT)
```

I removed the parameter `--phred64`, but it gives another error:

```
Saw ASCII character -126 but expected 33-based Phred qual.
terminate called after throwing an instance of 'int'
bowtie2-align died with signal 6 (ABRT)
```

Any idea?

Thanks

Answer: I struggled also several times with this kind of problem. Somehow (I don't know why) you may have a non-ascii character in your fastq file.

And for some reason (which I also don't know) bowtie2 explains that it has the ascii code "-126" which does not exist, as all ascii signs have to be positive.

It is relatively likely that you have non-ascii character in your quality lines of your fastq file.

You can search for a line containing a non-ascii character by using the following commands. If the file is gziped use this:

zcat FILE.fastq.gz | perl -e '$line = 1; while (<>) {if(/[^[:ascii:]]/) {print "LINE: $line\n$_";} $line++;}'

Otherwise this:

perl -e '$line = 1; while (<>) {if(/[^[:ascii:]]/) {print "LINE: $line\n$_";} $line++;}' -f FILE.fastq

Let's say you received such an output:

LINE: 21410096

CCCCCGGGGGGGGGG�GGGG

You non-ascii sign will be depicted as this question-mark symbol.

After receiving the responsible line you can "enlarge" the "surrounding area" by using sed:

sed -n '21410093,21410096p' FILE.fastq

So we got something like this as an output:

@sequence HEADER

CTGGCTGGGAAGGGGCTGGCT

+quality HEADER

CCCCCGGGGGGGGGG�GGGG

Such an corrupted entry can easily be removed by using sed. In this case the corrupted read was located at the lines
21410093 to 21410096, so the following command will work:

sed -n '21410093,21410096d' Fasta.fastq > repaired.fastq

This would be problematic if you have paired-end reads. Either remove the read entry at the same position in the corresponding fastq file, or replace the corrupting character in you initial fastq file. I would prefer the latter:

This can be done by using sed:

sed '21410096s/CCCCCGGGGGGGGGG�GGGG/CCCCCGGGGGGGGGGIGGGGG' FILE.fastq >FILE.repaired.fastq

(actually I'm surprised that this one worked).

I hope this helps.
Question: Hi,
I have a text file of 7500 genomic positions with their chr, start, end and want to get their nucleotide sequences. Can someone point me some tool or any thoughts on how to do it?

Thank you,

Answer: https://www.biostars.org/p/194822/#194987
Question: I am trying to extract only the columns I need from VCF preserving VCF structure, its header, its formatting. I am using bcftools. I tried doing:

```
bcftools annotate -c CHROM,POS,ID,REF,ALT,QUAL,FILTER,INFO/AF,INFO/AC,INFO/AN Holland.vcf -o Holland_selected_cols.vcf
```

But the output file just stays the same. Then I tried `query`:

```
bcftools query -f'[%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t%INFO/AF;%INFO/AC;%INFO/AN\n]' -H Holland.vcf -o Holland_selected_cols.vcf
```

But it does not preserve `VCF` header. What would be the right `bcftools` command for that?

Answer: Use `bcftools annotate -x` to remove all fields, except those you want to keep:

bcftools annotate -x ^INFO/AF,^INFO/AC,^INFO/AN,^FORMAT input.vcf
Question: Hello,

I know the following question have been asked many times by now, but the suggested solutions were not helpful.
I have a set of taxIDs and I want to convert them to taxonomies, but I do not want to download any extra databases. Is there any simple way to do this via terminal?

Answer: Using [Entrezdirect][1]:

$ esearch -db taxonomy -query "4932 [taxID]" | efetch -format nativemode xml | xtract -pattern Taxon -block "*/Taxon" -unless Rank -equals "no rank" -tab "\n" -element Rank,ScientificName
superkingdom	Eukaryota
kingdom	Fungi
subkingdom	Dikarya
phylum	Ascomycota
subphylum	Saccharomycotina
class	Saccharomycetes
order	Saccharomycetales
family	Saccharomycetaceae
genus	Saccharomyces

OR

$ esearch -db taxonomy -query "4932 [taxID]" | efetch -format native -mode xml | grep ScientificName | awk -F ">|<" 'BEGIN{ORS=", ";}{print $3;}'
Saccharomyces cerevisiae, cellular organisms, Eukaryota, Opisthokonta, Fungi, Dikarya, Ascomycota, saccharomyceta, Saccharomycotina, Saccharomycetes, Saccharomycetales, Saccharomycetaceae, Saccharomyces

[1]: http://bit.ly/entrez-direct