cannin/biostars_qa huggingface dataset:
Question: Hi I am trying to calculate the CAI (codon adaptation index) for a list of genes of interest.
I decided to use the RSCU values from the Codon Bias database (CBDB), for the relative organism (*B. subtilis* 168): http://homepages.luc.edu/~cputonti/cbdb/genera/bacillus.html#s30 .

I implemented a small python script to do so. While doing so, I couldn't decide between two methods tested them both: one is from Biopython SeqUtils (module CodonUsage, method cai_for_gene) and the other from the 'CAI' package (method CAI) (https://pypi.python.org/pypi/CAI). Since I was using the same reference RSCU values, and both claim to implement the method from "Sharp & Li, 1987, **NAR** (15)", I was expecting the same value for each gene from both packages. But this did not happen and it is clear at first sight as the results from the SeqUtils package are not included between 0 and 1 (in some cases are >1), which is in contrast with the definition of CAI, while for the 'CAI' package results were within the range.

I further investigated trying to regenerate the RSCU values and looking at both source codes.
The first thing I noticed is that from the same fasta file (i.e. same collection of genes, still downloaded from CBDB) I was getting two different RSCUs values on the two packages, with CAI-generated one matching the online posted values, while the values generated by SeqUtils were different. I assumed that the dictionary used by SeqUtils was not RSCUs values, thus I used the SeqUtils-generated RSCUs value to calculate the CAI for my genes, but still values were not allowed. Then looking at the source code I noticed that SeqUtils implementations differs a lot from the paper of Sharp and Li, unless it is a mathematical equivalence I do not know/understand. On the other hand the approach implemented in the CAI package is identical to the paper.

This given, I assume the CAI package is returning the correct values.

My question is thus: what is then SeqUtils calculating? I could not retrieve a lot of info about this package. Could there be an error? Why the two algorithms differ? And why the one from SeqUtils differs also from the paper? Anyone able to explain this?

P.S. I am of course assuming I used correctly both methods. Since I followed the relative guides, I feel reasonably safe in assuming so.

Answer: SeqUtils is actually not a correct implementation of the CAI metric (as of 6/22/18). I wrote the [CAI package](https://github.com/Benjamin-Lee/CodonAdaptationIndex) specifically out of frustration with SeqUtils.CodonUsage. In addition to being correct, the CAI package is faster and supports multiple genetic codes.

For more information, take a look at [the preprint][1] of the CAI software package paper.

Hope this helps and be sure to contact me if you run into any problems!

Benjamin Lee

[1]: https://github.com/openjournals/joss-papers/blob/joss.00769/joss.00769/10.21105.joss.00769.pdf "the preprint"
Question: I need to check which of the predicted de novo transcripts have polyA at 3'. Are polyA tails filtered at some step of tophat/cufflinks processing or this information is retained? If this info is retained, where do I find it. thank you

Answer: It's unclear what TopHat/Cufflinks has to do with the question as you are talking about de novo assembled transcripts. I would just start by using the grep command to get a feeling for whether there are any poly-A stretches. It is easy enough to write a script (Perl/Python/whatever) to look for them as well.

Often there are very very few reads with poly-A in HiSeq RNA-seq data, so you might not find much. For some reason there seem to be more poly-A tails left in MiSeq-produced RNA-seq data.
Question: I have a phased VCF file and I need to change the alleles of the SNPs so that the alleles coded as 0 represent the ancestral state while the alleles coded as 1 represent the derived state.

I got a list of ancestral and derive states for my SNPs in my VCF file and I want to change them using PLINK.

If I have a list of ancestral alleles in a file `ancestral_snps.txt` , if I use the following command will I get the wanted output?

`plink --vcf myvcf.vcf --a1-allele ancestral_snps.txt --recode vcf --out myvcfancestral.vcf`

Or is it `--allele-a2` ?

I assume that this change won't affect the phasing of my data, right? As I basically only recoded the alleles.

**EDIT:**

Following Christopher Chang's comment:

If my `ancestral_snps.txt` file looks like this, specifying which allele is ancestral:

C rs1001
G rs1002
A rs1003

If I use this command in `PLINK2`, will the ancestral alleles be coded as 0 in my phased VCF?

`./plink2 --vcf myvcf.vcf --ref-allele force 1 2 ancestral_snps.txt --recode vcf --out myvcfancestral.vcf`

Answer: - If you want to preserve phase information, you'll need plink 2.0; this is impossible in 1.9 due to the limitations of the core file format.

- "--ref-allele force ancestral_snps.txt" should work. Depending on the layout of the ancestral_snps file, you might need to specify some column numbers too; "plink2 --help ref-allele" provides usage details.
Question: Hi all,

I have read through many posts about insert size here. And see [a very good answer about the insert size][1].

It is still not so clear for me to understand insert size. I hope some experts can make it clearer.

As illustrated in [a good blog][2] and [a good anwser][3], the "insert size"=sequence between adapters (actually encompasses R1 and R2 as well as the unknown gap between them) and it is also known that [the ninth column of the SAM file (TLEN) represents the insert size][4]

However, here are some things I still don't understand.

First, in RNA seq data, if the alignments are spliced, and the TLEN reports the distance from the 5'-most to 3'-most position ([if my understanding is right][5]). So according to my understaning the TLEN number will include the possible introns which means the TLEN would be unsally longer than "actual insert size"?

Second, if we are mapping DNA sequences, then the fragment length and "insert size"/"template length" are the same?

Third, how Picard tools [CollectInsertSizeMetrics ][6] actually do to calculate the insert size distribution of a paired-end library, does it only use the TLEN or exclude possible introns?

Any answer to help me better ubderstand this conception will be greatly appreciated.

[1]: https://www.biostars.org/p/94246/
[2]: http://thegenomefactory.blogspot.de/2013/08/paired-end-read-confusion-library.html
[3]: https://www.biostars.org/p/106291/
[4]: https://www.biostars.org/p/94246/
[5]: http://samtools.github.io/hts-specs/SAMv1.pdf
[6]: https://broadinstitute.github.io/picard/command-line-overview.html

Answer: > So according to my understaning the TLEN number will include the possible introns which means the TLEN would be unsally longer than "actual insert size"?

Yes, the `TLEN` field won't always be terribly useful in RNAseq. When trying to compute the original fragment sizes it's best to not have spliced fragments. Back when we used tophat2 in our production pipeline, our "insert size estimation" step aligned to the transcriptome to avoid this problem.

> Second, if we are mapping DNA sequences, then the fragment length and "insert size"/"template length" are the same?

This will depend a bit on which fragment you're talking about. See the comment from https://www.biostars.org/u/18713/.

> how Picard tools CollectInsertSizeMetrics actually do to calculate the insert size distribution of a paired-end library, does it only use the TLEN or exclude possible introns?

I've never run that tool on RNAseq data, I'm not sure how useful it would be. I would expect that it's just summarizing the `TLEN` field, so I'd expect some absurdly high mean values.
Question: Hello,

I'm new here and and I need help please!

I need to download a list of all human genes with their respective gene symbol | chromosome | strand | transcription start site | Txen | and **Ensembl gene name**,

Actually I'm using UCSC table and I get something like this:

```
#hg38.knownGene.name hg38.knownGene.chrom hg38.knownGene.strand hg38.knownGene.txStart hg38.knownGene.txEnd hg38.kgXref.geneSymbol
uc001aaa.3 chr1 + 11873 14409 DDX11L1
```

And I want to know if it's possible to include ensembl gene symbol with UCSC table or with another method

Thank you in advance

Cherif

Answer: Follow these steps: http://kandurilab.org/bioinformatics/biostars/UCSC_tableBrowser_annotation.pdf
Question: I have a database of 10 distinct biological samples. Each of these samples was sequenced (RNA-Seq) using paired-end reads and 6 barcodes. Thus, for each of the ten biological samples, I have 12 .fastq files, with names like *ATACTC_1, *ATACTC_2, *GTGCTC_1, *GTGCTC_2...and so on.

I would like to follow this pipeline:

1. Analyze data quality with FastQC
2. Trim data with Trim Galore!
3. Align with TopHat2
4. ??? generate counts, differential expression analysis, etc.

Here is my question: at what point in this pipeline can I (should I) combine all of the .fastq data together? If each sample has 12 files associate it, at what point do I collapse the 12 files into a single file representing the RNA-Seq data for a single biological sample that I can analyze for counts in step #4?

I'm guessing I combine all of the .fastq files up front (re-multiplex?) with `cat file1...file12`. I could also do steps #1-2 or steps #1-3 completely, and then combine the output of step #3.

Thank you for any help you can provide! This board has already been tremendously helpful to me.

Answer: I think it's most typical to combine raw reads if they correspond to the same library and original sample (so merging technical replicates, with respect to the sequencer, and not biological replicates). [This paper][1] was one of the early ones to test and validate this assumption.

So for your outlined process, that could be anywhere in steps 1 to 3. Personally, I would make count tables for all the 6*10 runs separately and then combine the count tables at the very end, before any clustering or differential analysis (so in the middle of your step 4). This is because all the earlier steps can be performed in parallel, and you might save some time by processing many chunks at once.

And just for completeness, you'll need to combine the two matching paired end fastq files (e.g. `X_1` and `X_2`) before aligning. You might need to combine them or analyze them together for the adapter trimming, or possibly trim adapters for each read end separately.

[1]: http://www.ncbi.nlm.nih.gov/pubmed/18550803
Question: I've a table of blast results with ~5hits/query protein:

```
Protein Class
ProtA 1
ProtA 1
ProtA 1
ProtA 0
ProtA 1
ProtB 1
ProtB 1
ProtB 0
ProtB 0
ProtB 1
```

I would like to convert this into a feature vector matrix like this:

```
Protein Class1 Class2 Class3 Class4 Class5
ProtA 1 1 1 0 1
ProtB 1 1 0 0 1
```

Can someone suggest me an efficient way to do this, since I've ~2300k hits in the file.

Answer: Lets do it in Perl; Just change the filename in the first line, you should get the required output.

```
open FH,"makeintomatrix.txt";
my $dummy=; ##skip header
while()
{
chomp($_); # remove new line
#print "$_";
@array=split;
#print "$array[0]:$array[1]\n";
push @{ $HoA{$array[0]} }, "$array[1]"; ##hashes ofarray

}
print "Protein\tClass1\tClass2\tClass3\tClass4\tClass5\n";
foreach $k (sort keys %HoA) ##print the hash, with all the values
{
print "$k\t";

if("$#{$HoA{$k}}"==4)
{
foreach $tmp (@{ $HoA{$k} })
{
print "\t$tmp";
}
print "\n";
}
else
{
$fill=4-"$#{$HoA{$k}}";
#print "should be filled $fill times\n";
foreach (@{ $HoA{$k} })
{
print "\t$_";
}
for($i=1;$i<=$fill;$i++)
{
print "\t0";
}
print "\n";

}
}

close(FH);
```

`_____Data_______`

```
Protein Class
ProtA 1
ProtA 1
ProtA 1
ProtA 0
ProtA 1
ProtB 1
ProtB 1
ProtC 1
ProtC 1
ProtC 1
ProtD 0
ProtD 1
ProtD 1
ProtD 1
```

`_____Output_______`

```
perl makematrix.pl
Protein Class1 Class2 Class3 Class4 Class5
ProtA 1 1 1 0 1
ProtB 1 1 0 0 0
ProtC 1 1 1 0 0
ProtD 0 1 1 1 0
```
Question: Hey guys,

I am going to aligning reads from GSE35641 (http://www.ncbi.nlm.nih.gov/Traces/sra/?view=search_seq_name&exp=SRX119208&run=&m=search&s=seq), on ecoli reference genome (NC_000913). There are 6 sra files in GEO which the least size is more than 2G. Do you know any way to download the files faster? I need many hours to download a file of them in fasta format

Thank you biostars

Answer: Yep, install [Aspera][1] and run it as follows:

```
ascp -i $KEY/asperaweb_id_dsa.openssh -k 1 -QT -l 200m \
anonftp@ftp-trace.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR390/SRR390728/SRR390728.sra \ # url for the file
./ # save location
```

`$KEY` should point to the directory containing `asperaweb_id_dsa.openssh` (it will be in the installation directory). See instructions [here][2].

[1]: http://downloads.asperasoft.com/en/downloads/8?list
[2]: http://www.ncbi.nlm.nih.gov/books/NBK158899/#SRA_download.when_to_use_a_command_line
Question: I'm a newbie at snakemake.
I'm trying to implement the GATK FastqToSam in a rule.

I have the following and it works if I hard code the samplename into the shell part but I was wanting to get the samplename from the config file. I tried using the params option like so:

configfile: "config.yaml"

fc = config["flowcell"]
index = config["index"]
samplename = config["sn"]
lane=["L01","L02","L03","L04"]

rule all:
input:
expand(["result/gatk4/{fc}_{lane}_{index}_fastqtosam.bam"], fc = fc, lane = lane, index = index)

rule fastqToSam:
input:
trimmed1 = "result/trimming/{fc}_{lane}_{index}_1_500k_trimmed.fq.gz",
trimmed2 = "result/trimming/{fc}_{lane}_{index}_2_500k_trimmed.fq.gz"
output:
output="result/gatk4/{fc}_{lane}_{index}_fastqtosam.bam"
log:
"result/logs/fastqtosam/{fc}_{lane}_{index}.out"
benchmark:
"result/benchmarks/{fc}_{lane}_{index}.out"
container:
config["containers"]["gatk4"]
params:
RG="{fc}_{lane}",
sn="{samplename}"
threads: 8
shell:
"""
gatk --java-options "-Xmx8G" FastqToSam
FASTQ={input.trimmed1}
FASTQ2={input.trimmed2}
OUTPUT={output.output}
QUALITY_FORMAT=Standard
READ_GROUP_NAME={params.RG}
SAMPLE_NAME={params.sn}
LIBRARY_NAME=PCRFree
PLATFORM=MGI
SEQUENCING_CENTER=CG 2>{log}
"""

But I get the error:

WildcardError in line ** of snakefile:
Wildcards in params cannot be determined from output files. Note that you have to use a function to deactivate automatic wildcard expansion in params strings, e.g., `lambda wildcards: '{test}'`. Also see https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#non-file-parameters-for-rules:
'samplename'

Is there a way of getting the samplename from the config file and using it in the shell part of the snakemake rule?

Answer: This happens because you used the Snakemake notation for wildcards "{samplename}". Since you've already loaded the configuration file as a dictionary and stored the sample name value in the "samplename" variable, you can simply pass it as is:

```python
params:
RG="{fc}_{lane}",
sn=samplename # <------ HERE
threads: 8
shell:
"""
gatk --java-options '-Xmx8G' FastqToSam # <--- I also fix some quotation here the might bring you some problems.
FASTQ={input.trimmed1}
FASTQ2={input.trimmed2}
OUTPUT={output.output}
QUALITY_FORMAT=Standard
READ_GROUP_NAME={params.RG}
SAMPLE_NAME={params.sn}
LIBRARY_NAME=PCRFree
PLATFORM=MGI
SEQUENCING_CENTER=CG 2>{log}
"""
```
If I'm correct, this should work.

However, reading input parameters from the configuration file can be somewhat inflexible. I recommend following the Snakemake documentation and reading input information from an input table instead.
Question: Dear all,

I have downloaded the GRCh38 reference human genome in fasta format from the ftp.ncbi.nlm.nih.gov following the suggestions of this [post][1] and [this][2]. Now I need to generate a fusion genome thus I need the headers right. The headers of the GRCh38 are in this form:

>chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38
>chr2 AC:CM000664.2 gi:568336022 LN:242193529 rl:Chromosome M5:f98db672eb0993dcfdabafe2a882905c AS:GRCh38
>chr3 AC:CM000665.2 gi:568336021 LN:198295559 rl:Chromosome M5:76635a41ea913a405ded820447d067b0 AS:GRCh38
[...]
>chrUn_GL000218v1 AC:GL000218.1 gi:224183305 LN:161147 rl:unplaced M5:1d708b54644c26c7e01c2dad5426d38c AS:GRCh38
>chrEBV AC:AJ507799.2 gi:86261677 LN:171823 rl:decoy M5:6743bd63b3ff2b5b8985d8933c53290a SP:Human_herpesvirus_4 tp:circular

May I ask what is the format of the header? What the individual fields refer to?

Thank you

[1]: http://lh3.github.io/2017/11/13/which-human-reference-genome-to-use
[2]: https://www.biostars.org/p/318811/#334309

Answer: AC: Locus identifier / Accession number at Genbank & NCBI

gi: GI number (sometimes written in lower case, "gi") is simply a series of digits that are assigned consecutively to each sequence record processed by NCBI. The GI number bears no resemblance to the Accession number of the sequence record.

LN: Length of the contig

rl: role in the assembly, (chromosome, unplaced, EBV decoy)

M5: MD5 checksum (file integrity)

AS: Assemby version

SP: Species

tp: probably type, indicates if linear or circular

Anyway, to keep things compact I always delete everything after the first whitespace prior to building a genome.fa or index.
Question: Hi,
I have an RNAseq data from TCR alpha and bed chain single cell sequencing which I am trying to analyze. I looked at the post https://www.biostars.org/p/148371/ and also read the Mark Davis paper and the protocol, which I thought would be appropriate to follow. But I cannot either find the dataset from the Mark Davis paper not the tool they used (vdjfasta). I was writing to ask what is the latest tool that would be appropriate for assembling the reads and getting a summary the TCR clones from my dataset.
Thanks,
- Pankaj

Answer: Software for T/B-cell receptor repertoire analysis using non-amplicon samples:

1. [TraCeR](https://github.com/Teichlab/tracer) T-cell receptors, single-cell RNA-Seq
2. [VDjer](https://github.com/mozack/vdjer) B-cell receptors, RNA-Seq

P.S. Check out [T cell fate and clonality inference from single-cell transcriptomes (Nature Methods 2016)](http://www.nature.com/nmeth/journal/v13/n4/full/nmeth.3800.html) that describes your case.
Question: Hi All,
I am a newbie to sequencing technologies.

Is there a way to quantify the paired-end sequencing overlap. since completely overlapping the read pairs would be a waste of sequencing resources. A bit of overlap can be useful when doing the alignment but a small gap is optimal to maximize coverage.

Is there a way to check sequencing efficiency by using the alignment files. Because if there is lots of overlap we want to improve this.

Any help is appreciated.

Thank you.

Answer: Yes, you can do this with alignment BAM/SAM file. You can extract the record for concordant alignment `YT:Z:CP` from SAM/BAM file. Once you have concordant alignment, you can look for field 9 in SAM/BAM file. Field 9 (9th column) represents the fragment length of paired-end sequences.

From here, you can get the fragment length distribution. Based on your sequencing protocol, you should have the insert size for paired-end sequences. Then you should compare the fragment length with the insert size. If the fragment length is less than the sum of two reads, it means your paired sequences are overlapped. Here, you can plot the histogram of fragment length distribution.

Note: concordant alignment record will give you the alignments which are within the given insert size
Question: I have a fasta file of yeast orfs (from the [SGD][1] database) and I need to find out (predict) secondary structure/classification for all of them (so ~6K queries and doing this by hand on numerous available websites is not an option).

Is [CATH][2] database a place that can give me this info? (so far I struggle to recover anything that can help me from the website but perhaps I am missing something), if not what's the best way/software to go about it? perhaps a piece of code in R or Python.

Surely, these predictions have been done before and there is a database somewhere? -- yeast is too much of an important organism, with a lot of information about its genome.

[1]: http://www.yeastgenome.org/
[2]: http://www.cathdb.info/

Answer: Okay, so after some searching around I have settled on [s2d][1] software, which might be a bit fiddly to set up but after the initial pain is very easy to use, with nice parseable output.

I've also been recommended [Jpred4][2] -- this one's easy to set up and use.

Both programs offer batch submissions and with s2d I was done with my ~6K proteins in a couple of days (but with some planning I could have done it even faster).

[1]: http://www-mvsoftware.ch.cam.ac.uk/
[2]: http://www.compbio.dundee.ac.uk/jpred/
Question: Hi all, I am new to the BioInformatics, and quite a beginner in programming languages. Can anyone suggest me some sources where I can at least learn 50% of the scRNA seq data analysis? I am familiar with C language, and I know a little bit of molecular biology too.

Answer: Yes, Seurat would be one of the starting points. Beside the tutorials offered on Seurat web site, a while I have posted some R code on Seurat github page : https://github.com/satijalab/seurat/issues/1193 (hope it is helpful)
Question: In the [Adapterama paper][1], they are talking about a way to tag individuals within a library preparation for sequencing.

At some point, they explain their primer and show this:

` iTru5_01_A: 5';AATGATACGGCGACCACCGAGATCTACAC tag;ACCGACAA 3';ACACTCTTTCCCTA*C`

What is the star for in this sequence:

`ACACTCTTTCCCTA*C`

This appears [in other protocols][2].

[1]: http://www.biorxiv.org/content/biorxiv/early/2016/06/15/049114.full.pdf
[2]: https://github.com/eacooper400/GenDivergence/blob/master/dRAD_LibraryPrep_Protocol.docx

Answer: It's some sort of modified nucleotide: phosphorothioate bond

See this page from [IDT][1]

[1]: https://www.idtdna.com/pages/decoded/decoded-articles/core-concepts/decoded/2011/09/12/phosphorothioates
Question: Hi,

I have been doing some ancient DNA bioinformatics for a while now and have been using the bwa aln command.

When using this command with ancient DNA it is recommended to disable seeding. This slows down alignments somewhat but allows the damaged DNA sequences to be aligned, the characteristic damage is a C->T transition at the 3' end of the sequenced fragments. For instance I use the parameters `-l 1024` (to disable seeding) `-o 2` (increase gap open slightly) and `-n 0.03` (reduce the edit distance slightly). These parameters were taken from this paper

http://www.biomedcentral.com/1471-2164/13/178

**My question:**

Are there any specific settings required to use bwa mem with ancient DNA? Does bwa mem use seeding at all?

Thanks in advance

Answer: I had really good experiences with `bwa mem` when using it on dataset with high error rates. It far outperformed other tools. See this:

https://www.biostars.org/p/117225/#117259

Your problem seems to be very similar in nature in that it introduces mismatches that are similar to sequencing errors. From what I understood from the paper `bwa mem` was designed to be robust to this type of errors.
Question: Hi everybody,

I have proteome data for a non model organism, and I would like to perform gene onthology (GO) enrichment analysis.

I would like to perform the enrichment for the up-regualted proteins using the all proteins encoded in the genome as background dataset. I downloaded the GO annotation for each protein directly from Uniprot and now I am kind of lost among the many tools available for doing the analysis.

I am looking for web tools or R-Bioconductor tools more appropriate form my purpose, because it seemed to my that most of the tools deal with model organism.

Any advice will be greatly appreciated.

Answer: An update.

I used the Bioconductor package [topGO][1] and I think it is a great tool for performing enrichment for non model organisms.

[1]: http://www.bioconductor.org/packages/release/bioc/html/topGO.html