cannin/biostars_qa huggingface dataset:
Question: I have two files:

1. An annotation file in gff format
2. A fasta formatted sequnce file for a whole genome.

I want to have a resulting file : myGenomeCDS.fasta

How can I go about doing this via python or Tools ?

Thanks

Answer: If your gtf file has the annotation of CDS's you can extract CDS sequences using `gffread` in this way:

gffread -g genome.fa -x CDS.fa annotation.gtf
Question: I have some RNA sequencing reads to align to the human reference genome. I found the genome FASTA files on both [GENCODE][1] and [ENSEMBL][2]: `GRCh38.p13.genome.fa.gz` and `Homo_sapiens.GRCh38.dna.toplevel.fa.gz`
But after unzipping them, I found that they are 3.1G and 60G respectively. Why is that? And which one should I use? (considering the purpose of the project is to detect gene fusion from the sequencing reads).

[1]: http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/
[2]: http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

Answer: https://www.biostars.org/p/487975/
Question: Hi there,

I was working for 5 years with NGS data and CNV detection, but now I have to perform CNV calling in .CEL files from Affymetrix v6 arrays. At first, I have no idea how to operate with such files (should I download some converter to vcf or txt or other format?). At second, I need to perform tumor CNA calling there (of course I have tumor-normal data). Which tool is the best, taking into account purity and clonality of samples?

Thanks a lot for any hints!

(of course I googled for papers, but I am a complete newbie in this field and a bit overwhelmed, everybody says 'oh we are the best buy our software', but somehow I dont trust most of them and also I work in academia and looking for something which is free for research)

Answer: Around the time when NGS was becoming very popular, the then Affymetrix CEO, Frank Witney (I managed to attend a meeting with him once), still firmly wanted to focus the company on microarrays. At the time, they [Affymetrix] had been pushing their OncoScan product. Prior to this was the Affymetrix SNP 6.0 array (SNP 6.0 array), which has probes for both SNPs and CNVs, as you know. The leading method that emerged for calling copy number segments in this SNP 6.0 array was Circular Binary Segmentation (CBS). The Aroma project did a lot of work in producing materials and guides for, for example, calling somatic copy number segments between tumour and normal. Take a look: https://aroma-project.org/vignettes/

Their code should still work. You will require the CDF for the SNP 6.0 Array, though. If you do not have it, it should be available on ThermoFisher's website (they bought Affymetrix). I may also have it on my disk somewhere.

Kevin
Question: Hi everyone,

I am new to the NGS. For training, I used paired-end data from public files, and usually, there were only two fastq files for the two reads.

Now, I receive the fastq files for my ChIP-seq studies. We did 50 bp paired-end. And for one PE file, there are four fastq files. For example:

R1F=IN1_S1_L001_R1_001.fastq
R1R=IN1_S1_L001_R2_001.fastq
R2F=IN1_S17_L002_R1_001.fastq
R2R=IN1_S17_L002_R2_001.fastq

I am now a bit confused about how to handle these 4 data.

I tried to do bowtie2 paired-end mapping for R1F as -1 and R1R as -2 but it seems not working.

Bowtie2 notification:
Note that if files are specified using -1/-2, a file cannot
also be specified. Please run bowtie separately for mates and singles.
Error: Encountered internal Bowtie 2 exception (#1)

I am also unsure about the use of Paired-end merge tools.

Please help. Thank you

Answer: Hi! I think what you have there is 2 fastq files for the same sample, each one coming from a different sequencing lane (note that some of the files have L001 and some others L002). What you need to do is to concatenate the sequencing lanes into one file and then align with the resulting file.

cat IN1_S1_L001_R1_001.fastq IN1_S1_L002_R1_001.fastq > IN1_S1_concatenated_R1_001.fastq
cat IN1_S1_L001_R2_001.fastq IN1_S1_L002_R2_001.fastq > IN1_S1_concatenated_R2_001.fastq

The resulting files are the 2 mates (R1 and R2) that you want to feed into your favorite aligner.

Hope it helps!
Question: Hi,

I have these lists of genes

> dim(b)
[1] 4866 2
>
How I can extract the small list (bb) from big on by the same order of genes???? I mean the order of genes from extracted file should be the same with small file(bb).

Answer: b and bb in OP doesn't share any gene. What is being asked, is not possible when values do not exist (between two data sets). Look at the following example and see if this is what is required:

> b=data.frame(genes=paste("gene", sample(10), sep="_"), expn=round(rnorm(10,1,4),2))
> bb=data.frame(genes=paste("gene",sample(5), sep="_"))
> b
genes expn
1 gene_6 -0.78
2 gene_10 4.65
3 gene_9 1.86
4 gene_3 2.39
5 gene_1 0.34
6 gene_2 2.01
7 gene_8 -4.51
8 gene_7 -7.61
9 gene_5 1.50
10 gene_4 3.12

> bb
genes
1 gene_2
2 gene_1
3 gene_3
4 gene_4
5 gene_5

> b[match(bb$genes,b$genes),]
genes expn
6 gene_2 2.01
5 gene_1 0.34
4 gene_3 2.39
10 gene_4 3.12
9 gene_5 1.50
Question: I'm trying to edit an MSA (Multiple Sequence Alignment) file generated by ClustalW, to trim sequences before the consensus one, using BioPython.
*xxx refers to other bases not relevant here*

Here's the example I/O :

INPUT

ITS_primer_fw --------------------------------CGCGTCCACTMTCCAGTT
RBL67ITS_full_sequence CCACCCCAACAAGGGCGGCCACGCGGTCCGCTCGCGTCCACTCTCCAGTTxxxxxxxxxxxxxxxx
PRL2010 ACACCCCCGAAAGGGCGTCC------CCTGCTCGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxx
BBF32_3 ACACACCCACAAGGGCGAGCAGGCG----GCTCGCGTCCACTATCCAGTTxxxxxxxxxxxxxx
BBFCG32 CAACACCACACCGGGCGAGCGGG-------CTCGCGTCCACTGTCGAGTTxxxxxxxxxxxxxxxx

EXPECTED OUTPUT

ITS_primer_fw CGCGTCCACTMTCCAGTT
RBL67ITS_full_sequence CGCGTCCACTCTCCAGTTxxxxxxxxxxxxxxxxxxxx
PRL2010 CGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxxxxxxx
BBF32_3 CGCGTCCACTATCCAGTTxxxxxxxxxxxxxxxxxxx
BBFCG32 CGCGTCCACTGTCGAGTTxxxxxxxxxxxxxxxxxxxx

The documented code for `AlignIO` describes just a way to extract sequences by treating the alignment as an *array*. In this example

align = AlignIO.read(input_file, "clustal")
sub_alignment = align[:,20:]
I was able to extract a subalignment made by all the sequences (:) starting from the 20th nucleotide. I'm looking for a way to replace the `20` in the example with the position of the first nucleotide of the consensus sequence.

Any answers including some cline software to trim easly as requested are well accepted. Will be great if python coded or for UNIX.

Answer: This should work. Simply finds the index of the first column that doesn't contain a gap character.

import sys
from Bio import AlignIO

# Usage: python scriptname.py alignment format

aln = AlignIO.read(sys.argv[1], sys.argv[2])

for col in range(aln.get_alignment_length()):
if not "-" in aln[:, col]:
position = col
print("First full column is {}".format(col))
break

print("New alignment:")
print(aln[:, position:])
Question: I have a wrapper python script that used RPy2 and I am calling native R functions and script, but I want to run R in a multithreaded environment. At the moment I am able to run R in multi process mode. Any ideas how this can be done?

Answer: There are multiple ways to do this, but one of the most common is to use packages like 'parallel', 'doMc' and 'foreach' to split out specific loops onto multiple cores in a map/reduce like manner.
Question: Hi,
Just like we have FastQC for QC of WGS and exome-data, is there a tool for QC of targeted-gene panel studies? The reason I ask is because some of the metrics from FastQC such as GC content may not hold for a small region captured in a focused NGS panel.

Thanks!

Answer: Check out CollectHsMetrics in Picard

https://broadinstitute.github.io/picard/picard-metric-definitions.html#HsMetrics
Question: Hi all,

Can you please explain to me the relation between Wald's test and Negative binomial generalized linear model? As for my understanding, the count data is modeled using negative binomial generalized linear model after which Wald's test is applied to figure out whether a particular gene is significant or not. Please correct me if I'm wrong.

Answer: RNA-seq raw count data 'naturally' follows a negative binomial distribution (Poisson-like), so, the DESeq2 authors model the data as such. By 'model the data', we merely imply that we build a regression model of the data such that we can make statistical inferences from it [the data].

So, after normalising the raw counts, the following occurs:

For each gene, a logistic regression model with the negative binomial as family is fit:

require(MASS)
gene1.model <- glm.nb(gene1 ~ CaseControl + ..., data=MyData)
gene2.model <- glm.nb(gene2 ~ CaseControl + ..., data=MyData)
*et cetera*

Once we have modeled each gene, a simple way to derive a P value for each model coefficient (i.e. CaseControl, etc) is by applying the Wald Test and selecting the coefficient of interest:

require(aod)
wald.test(b=coef(gene1.model), Sigma=vcov(gene1.model), Terms=c(2)) #term '2' would be CaseControl

The Wald test is a standard way to extract a P value from a regression fit.

Kevin

NB - this is ***not*** the exact code used by DESeq2, of course. This is just giving you a broad overview with some simple R functions. For one, DESeq2 models dispersion in addition to everything that I have mentioned above, and the Wald test is not used in each case to derive p-values in DESeq2.
Question: Hi folks:

I am trying to get raw RNASeq data from ENA and wonder how to know the length of the sequences in those FASTQ files without downloading them?

In the query design I follow their [guide][1], and use the "`nominal_length>XXX"` in the query, but this filtering mechanism doesn't work...

Thanks!

[1]: http://www.ebi.ac.uk/ena/browse/search-rest

Answer: The Paired nominal length does not seem to work in the ENA query builder [here][1] or at least it doesn't autofill into the query box, but you can add it to the URL string. For example, in the query builder, select domain = read and Library layout= Paired and optionally a taxon name like Bacteria , then you get this query.

tax_eq(2) AND library_layout="PAIRED"

Hit search and then add nominal length to the URL and get this [URL string][2] and query

tax_eq(2) AND library_layout="PAIRED" AND nominal_length>300

The nominal_length is just a filter, so you may need to check read experiment XML for actual length and stdev.

UPDATE : add `&result=read_experiment&display=xml&download=txt` to the URL and parse `NOMINAL_LENGTH` from [the summary XML][3]

[1]: http://www.ebi.ac.uk/ena/data/warehouse/search
[2]: http://www.ebi.ac.uk/ena/data/warehouse/search?query=%22tax_eq%282%29%20AND%20library_layout=%22PAIRED%22%20AND%20nominal_length%3E300%22&domain=read
[3]: http://www.ebi.ac.uk/ena/data/warehouse/search?query=%22tax_eq%282%29%20AND%20library_layout=%22PAIRED%22%20AND%20nominal_length%3E300%22&domain=read&result=read_experiment&display=xml&download=txt
Question: Hello,

I am trying to generate a .vcf file and call variants in RNA-Seq data, as per the GATK3 tutorial: https://software.broadinstitute.org/gatk/documentation/article.php?id=3891

I can get all the way through making a filtered .vcf file, but the resulting file has no entries - i.e., the .vcf file is only 47 lines, and all of that is header information. The columns are empty. The unfiltered .vcf file in the HaplotypeCaller is also empty.

My questions are:

- First - as a sanity check - this is an unlikely outcome, right? It's not that my RNA-Seq coverage was too shallow (~30,000,000 paired end reads) and therefore we are unable to detect any variants?
- Second - if we should, indeed, see some variants in that file - since it's hard to tell if there are variants at all until you actually generate a .vcf file (i.e., you can't look in a bam file and decide that variants exist) - what parameters should I tweak in order to try and get variants in my final .vcf file? I'm not sure what intermediate steps I can look at to troubleshoot where in the pipeline things are going wrong.

I didn't perform Indel Realignment (step #4) or base recalibration (step #5) because it was unclear to me from the tutorial if this was necessary for RNA-Seq. Right before step 6, I also had to use ReorderSam and BuildBamIndex in order to make the intermediate bam files generated match the ordering of the .fa file I am using as an index. I'll work on getting base recalibration working next, but if anyone has additional suggestions or answers to the above, I'd appreciate it.

Thank you for your help!

Kristin

Answer: Thank you everyone for providing valuable feedback that helped me narrow down the world of what could possibly be going wrong and identifying the key issue.

Ultimately, the problem was that I hadn't noticed that at some point the second alignment had failed, and I was trying to run all these commands on an empty .bam file. Thus - the problem wasn't the .dict at all - but the files that were being used for calling.

I'm a little confused about how I missed this because I had been using samtools to take a peek at every file and make sure they looked fine - but this troubleshooting process has taken so many iterations that I must have missed something.

I was inspired by this post, where someone encounters a similar issue: https://gatkforums.broadinstitute.org/gatk/discussion/10130/gatk-3-8-dictionary-cannot-have-size-zero

Now that I went back and make sure the first and second alignment steps produced .bam files, I do not get the .dict error, and my .vcf files have variants written in them.

Thank you all!
Question: Is there a good way to distinguish ./. from 0/0 in a merged vcf? For example, a tool that goes back to the bam files for missing genotypes and checks if it's homozygous reference or a NO CALL? I would imagine this would be important to do, especially in 30x WGS where there are a lot of low coverage areas that lead to no calls.

Although, I guess if you have the bam files maybe the best thing to do is joint variant calling?

Answer: I just wrote http://lindenb.github.io/jvarkit/FixVcfMissingGenotypes.html It takes a merged VCF and uses the original Bams to fill the missing genotypes. . If the number of reads is greater than **min.depth**, then the missing genotype is said hom-ref. This tool is very new/alpha, use the github issue tracker to tell me if there is a problem please.

**Usage:**

```
$ yourtool-mergingvcf 1.vcf 2.vcf 3.vcf > merged.vcf
$ find ./ -name "*.bam" > bams.txt
$ java -jar dist/fixvcfmissinggenotypes.jar -f bams.txt merged.vcf > out.vcf
```
Question: I'm looking for an API where I can enter the common name of a plant or organism into it, and it returns the scientific/latin name and vice versa. I see that ITIS has an API, they have a list of examples of how to use it here:
https://www.itis.gov/ws_searchApiDescription.html

They have a search by common name feature, but it doesn't actually return the latin name like it says its supposed to, heres what it returns:

https://www.itis.gov/ITISWebService/services/ITISService/searchByCommonName?srchKey=american%20bullfrog

Morning Glory Pellicia
English
706771

They have an API for getting common and latin names from the TSN:
https://www.itis.gov/ws_tsnApiDescription.html
but I don't know the TSN for the plants. Is there a way to use this API to convert between common and latin names, or is there another API or tool that does it?

Answer: With [NCBI Entrez Direct][1]:

esearch -db taxonomy -query "Morning Glory Pellicia" \
| efetch -format docsum \
| xtract -element ScientificName
Pellicia dimidiata

You can also do it with [eutils][2].

[1]: ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/
[2]: https://www.ncbi.nlm.nih.gov/books/NBK25501/
Question: where to get the file of hg19 exon, intron, UTR region ?

I read lots of post in biostar. However it's been a long time post.

Answer: Be aware that UTRs are part of exons. So when you get UTR regions from a GTF/GFF, they will overlap with the regions annotated as exons in the GTF. Perhaps what you want is not things annotated as exons, but things annotated as CDS?
Question: Hi, I used blasr aligner to align a pacbio dataset against a reference. I used *.bax.h5 files for the alignemnt. Since there were three *.bax.h5 file for a single sample I used the command as follows:

blasr m160128_145513_42215_c100901562550000001823208404291605_s1_p0.1.bax.h5 m160128_145513_42215_c100901562550000001823208404291605_s1_p0.1.bax.h5 m160128_145513_42215_c100901562550000001823208404291605_s1_p0.1.bax.h5 query.fasta -sam -unaligned unaligned.reads -nproc 8 -out output.sam

I converted the sam to bam and sorted it. Now, when I check the bam file for the counts of read mapped/unmapped, I get something like this:

samtools view -c -F 4 output.sorted.bam

12978

but when I use a different approach of counting reads I get 3 different values:

samtools view -F 4 output.sorted.bam | wc

12978 298494 9212184

I am not able to understand if these 3 values belong to individual *.bax.h5 I originally used in the mapping or this something else?

Or

is there any other way to use 3 bax.h5 file together in the same mapping step?

Thanks

Answer: `man wc`

12978 is the number of lines and exactly what `samtools view -c ...` reported.

298494 is the number of words. This is useless.

9212184 is the number of characters. This is useless.

You probably wanted `wc -l` rather than just `wc`.
Question: Hi,

I don't know why I can't sort my bam files

dhcp179185:unknown$ samtools sort SRR1593424.bam SRR1593424_sorted
[bam_sort] Use -T PREFIX / -o FILE to specify temporary and final output files
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-n Sort by read name
-t TAG Sort by value of TAG. Uses position as secondary index (or read name if -n is set)

Answer: Do

samtools sort -o SRR1593424_sorted SRR1593424.bam
Question: I have a fastq files and I found many reads with only N. I want to filter only those reads which contains N.
Any one liner for that?

Answer: To separate reads with Ns using BBMap:

bbduk.sh in=reads.fq out=readsWithoutNs.fq outm=readsWithNs.fq maxns=0

If you have, say, 100bp reads and only want to separate reads containing all 100 Ns, change that to "maxns=99".
Question: Hi guys,
I have a folder with around 300 `.bam` files. Each `.bam` file is a lane of a sample and hence 4 lanes make a sample. I would like to merge the .bam files of the four lanes in a single one by `_*S*_` where S is followed by a number that represent the number of the sample (e.g. `my_experimet_xxx__L001_S1_stimulated_Aligned.bam`).
Suppose I have 75 samples, i.e. `${1...75}`.

Can anyone help me please?

The line I use to merge normally is the following:

samtools merge S1_merged.bam *bam

Thank you in advance

Answer: I think this should work for your issue:

create a file name samtools_merge.sh

$mkdir merged

for L1 in *_L001_*.bam
do
echo $L1
L2=`echo $L1 | sed 's/_L001_/_L002_/'`
L3=`echo $L1 | sed 's/_L001_/_L003_/'`
L4=`echo $L1 | sed 's/_L001_/_L004_/'`
merged=`echo $L1 | sed 's/_L001_/_merged_/'`
samtools merge ./merged/${merged} ${L1} ${L2} ${L3} ${L4}
done

This will iterate over each unique sample with L001 somewhere in the name and store other variables by replacing L001 with L002,003,004, and do this for each sample. Then it will run samtools merge on all 4 lanes, then do the same for the next sample until it has gone through all the samples. The filename output will be the same as the sample name, but will substitute lane information with "merged".

Just run this in your directory with all the bams and you should have merged bams in the dir "merged".
Question: Hello everyone

I have the following loop to get the TaxID from a list (aaa.txt) of Assembly Accession.

head aaa.txt
GCF_016924235.1
GCF_003999975.1
GCF_002993365.1
GCF_017084525.1
GCF_006496635.1
GCF_017104785.1

The loop:

cat aaa.txt | while read p; do
echo $p;
esearch -db assembly -query $p | esummary | xtract -pattern DocumentSummary -def "NA" -element AssemblyAccession,Taxid > out.txt;
done;

I can't figure out why the loop stops after the first line iteration.

Thank you

A

Answer: As @cpad0112 suggested, adding `< /dev/null` to your `esearch` is the way to go if you want to use the while loop.
I recommend skipping the while loop altogether and use `epost` as shown below. It is quicker.
```
$ time cat accs.txt | while read p ; do esearch -db assembly -query $p < /dev/null | esummary | xtract -pattern DocumentSummary -def "NA" -element AssemblyAccession,Taxid ; done > out.txt
real 0m5.361s
user 0m1.828s
sys 0m0.553s

$ time epost -db assembly -input accs.txt | esummary | xtract -pattern DocumentSummary -def "NA" -element AssemblyAccession,Taxid > out.txt
real 0m1.567s
user 0m0.375s
sys 0m0.150s
```
If you have a lot of accessions and you are concerned about hitting the rate limit, I suggest you create an NCBI API key. See https://support.nlm.nih.gov/knowledgebase/article/KA-05316/en-us and the "Programmatic Access" section of https://www.ncbi.nlm.nih.gov/books/NBK179288/ for more details.

Finally, if you don't want to deal with EntrezDirect at all because it still takes a long time to process thousands of accessions you have, you can obtain this information from the assembly report files located on the FTP here: ftp://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS. For example, the file `assembly_summary_refseq.txt` is a tab-delimited table with 23 fields including assembly accession and taxid.

You should be able to join on the assembly accession to this table:
```
$ join -j 1 <(sort accs.txt) <(grep -v '^#' assembly_summary_refseq.txt | sort -k1,1 -t $'\t') -t $'\t' | cut -f1,6
GCF_002993365.1 2079529
GCF_003999975.1 2486853
GCF_006496635.1 2027405
GCF_016924235.1 2811233
GCF_017084525.1 2812560
GCF_017104785.1 2703894
```
Note, `assembly_summary_refseq.txt` contains only the latest assembly accessions. Older data are in `assembly_summary_refseq_historical.txt` file located in the same FTP path. If the input list has a mix of current and old assemblies, the best way to go is to first concatenate the two files and then do the `join`.
Question: Is anyone aware of a tool which allows one to apply the concept of "mutation hotspots" to 3D protein structures? It is common to define mutation hot spots at the DNA level where we define the hotspot as a region of DNA sequence with recurrent mutations. But what if we want to know if some mutations are "geographically" co-located in the 3D protein structure (e.g., in the same pocket) but not necessarily "sequentially" co-located. Even without pre-defined concepts of "pockets" if a protein structure is known then it should be possible to define some concept of physical distance between residues. Any suggestions?

Answer: It looks like you can also visualize mutations for certain genes in a 3D context using [CBioPortal][1]. Enter a gene such as ESR1. The select the 'Mutations' tab. Then hit the '3D Structure' button beside the mutation diagram (loliplot) that shows the 1D location of observed mutations along the amino acid sequence. This opens a side panel that shows a 3D model for each PDB chain that is available for your protein of interest. You are also give some ability modify the 3D visualization and the way mutations are annotated on it. Very slick interface.

[1]: http://www.cbioportal.org/public-portal/
Question: Hello.
In tool 'AddOrReplaceGroups''s arguments - --RGPL, --RGLB, --RGPU, --RGSM,
what do they mean by each one?

In case of. RGBL, should i set a name of read group library now or I should fill in exsisting name?

And I wanna know exactly what other things mean.
Please explain in detail.
Thank you in advance.

Answer: See the following for more guidance
https://software.broadinstitute.org/gatk/documentation/article.php?id=6472