cannin/biostars_qa huggingface dataset:
Question: hello, I am trying to run gatk cnv caller tutorial and I am facing an error with respect to reference fasta file. I have installed gatk using conda.
The code I am trying to run is:

gatk PreprocessIntervals \
-R ~/tutorial/GRCh38/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
--padding 0 \
-L ~/tutorial/tutorial_11684/gcnv-chr20XY-contig.list \
-imr OVERLAPPING_ONLY \
-O ~/output/chr20XY.interval_list

and the error i am getting is:

> A USER ERROR has occurred: The specified fasta file (file:///home/username/tutorial/GRCh38/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta) does not exist.

Please give some suggestions if there is any error in the code.

Answer: Have you checked to see if the file is readable? You can use

ls -lh
If it isn't readable, you can make it so using chmod. For example, you could use

chmod a+r
Question: Dear Biostars folk,

I have been trying to come up with an efficient way to return all the NON-hits from a blastn query. Basically, I need to blast all reads in an Illumina library against one another, and only find the reads in set A that have not hit in B. Although I have a solution, it seems to be really slow, and I wonder if you could help me out. Here's what I did so far:

makeblastdb -in ${B} -dbtype nucl
blastn -query ${A} -db ${B} > blast_results.dat
cat blast_results.dat \
| grep 'No hits' -B 6 \
| grep 'Query' \
| cut -d '=' -f2 \
| cut -d ' ' -f2 \
> unique_read_ids.dat
cat unique_read_ids.dat \
| while read line; do sed -n "/${line}/,/>/p" $query \
| sed '$d'; done \
> result.fn

Now... this totally works, but it takes about a week to finish, although the blastn-part only take 1-2 days. So that last bit, where I fish out all the non-hits, it probably not so efficient.

Is there any way to let blastn return NON-hits, instead of all queries / all hits? If so, how does that work? I've tried setting something up with the e-values, but it always seems to be a minimal value...

Thanks for your time!

Answer: For anyone returning to this in the future, there's a biopython tutorial on doing exactly this, in a much more efficient way: https://biopython.org/wiki/Retrieve_nonmatching_blast_queries
Question: Hi,

I have two bam files from ICGC. During mutation calling, vardict/mutect/freebayes/varscan split bam to by its chromosomes. Those GL00 parts also appear just like chr1-2-3-4-5-6-7 appear. I check what they are and found they are contigs?!?(forgive my ignorance) How should I treat them ? Will they cause a problem during variant calling output interpretation? I read one previous thread about this but there were not any information about how to treat them. I need your help.

Best,

Tunc.

Answer: You can check http://genome.ucsc.edu/cgi-bin/hgGateway where they provide details for every assembly.

This was addressed on the UCSC genome support forum:
http://redmine.soe.ucsc.edu/forum/index.php?t=msg&goto=8701&S=5061c99adf3f4c4ee90f0ea362d838c7

> chr*_random - are called "unlocalized" sequences. The chromosome is
> known, but not the location on the chromosome. chrUn_* - are called
> "unplaced" sequences. They probably belong to the sequenced genome,
> but placement is unknown at this time.
> The GL numbers (or other types of numbers) in these names are the
> genbank identification numbers which can be used in a nucleotide
> search at Entrez. For example: chr1_GL456211_random - unlocalized
> sequence belonging to chr1, NCBI identification: GL456211
> chrUn_JH584304 - unplaced sequence, NCBI identification JH584304
Question: Dear Biostar members,

I am involved in analyzing data from NGS technologies. As you may know, raw data files are in fastq format. And the fact is that these files are huge (several gigabytes each).

As a first step, we naturally want to make some quality controls of the sequencing step. My team and I are aware of some tools like Picard or bamtools that give some information and statistics, but they do not allow to extract all the information we want. So, at some point, I am afraid we will have to parse fastq files ourselves.

Which programming language would you choose to do the job (ie : read and process a 10GB file) ? Obviously, the main goal is to make the processing as fast as possible.

awk/sed

C

Perl

Java

Python

Do you think that we will be I/O bound anyway, so using a language or another makes no big difference ?

Would you read the file line by line or load chunks of hundreds of lines (for instance) ?

Many thanks.
T.

Answer: The short answer is: it depends a lot on what exactly it is you need to do with the contents of the files.

Regarding being I/O bound, I think it is quite unlikely unless you have the data on very slow disks. Even a single 5400rpm disk can read in the order of 50MB/s, which means that reading a 10GB file will not take much over 3min. Since you are working with large NGS datasets, I would assume that you have a RAID system, in which case the combined speed of multiple disks will almost certainly exceed how fast you can process the data.

When it comes to speed, nothing beats a good C or C++ implementation. However, there is a huge price to pay in terms of development time. If what you do is pretty basic string operations, both Perl and Python should do fine, and the savings in development time are most likely well worth the slightly longer run time compared to a C or C++ implementation.

In my experience, Perl has a more efficient implementation of regular expressions. Conversely, if you need to load all the data into memory, the Python data structures appear to be more memory efficient.

Awk is the only of the languages that you list that I would advise you against. It is very convenient for writing short command-line hacks, but speed is not one of its virtues. I don't have any experience with using Java for bioinformatics tasks.
Question: hi everyone,

I have just started learning genomics as a part of my bioinformatics degree and I've been introduced to linux for handling fastq files and using fast qc. can anybody suggest some good learning resources which is more inclined towards linux for genomics.

Answer: [The Linux Documentation Project][1] has many resources to read

[1]: https://tldp.org/index.html
Question: I have fasta file namely `119XCA.fasta` as shown below,

>cellulase
ATGCTA
>gyrase
TGATGCT
>16s
TAGTATG
I need to remove all the fasta headers, keep the sequences one by one and need to write file name as a fasta header. The expected outcome is shown below,

>119XCA
ATGCTA
TGATGCT
TAGTATG
I have used the following script `sed '/^>/d' foo.fa > out.fa` which remove the fasta headers but, i do not know how to manage to write file name as a header. Therefore, please help me to do the same.

Answer: Assuming you're using BASH, use `basename` to get the filename with no PATH. Like:

filename=$(basename -i file | cut -d'.' -f1)

Then you could replace it using sed

sed -i "s/^\>.*$/$filename/" your.fasta

Remember to use double quotes to use variables in sed.
Question: Hi,

I am using the gencode v19 gtf file for annotation of my chipseq data ( source - ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/) and the description of the format of the file (http://www.gencodegenes.org/gencodeformat.html) mentions that there is an optional field called transcript_support_level. However I cannot find that information for any of the records in the gtf file nor does any other data file in the release has that information.

Surprisingly though, in the UCSC browser, for all gencode transcripts, the transcript support level has been mentioned ( which to me seems like that there should be a file somewhere that contains this tag). The gtf files from UCSC gencode track/Ensembl FTP also do not contain this record.

Has anybody used this field before? and if so how? Does anybody have an idea on where can I source the data sheet that would contain the transcript support level information?

Thanks!

Answer: We brought them in for the most recent Ensembl/GENCODE release, which is Ensembl 77/GENCODE 21, so they will not be available in Ensembl 75/GENCODE 19. Also, at present they are not being included in GTF files.
Question: Hi,

We have run a pilot RNA-Seq study and I used `edgeR` package to obtain differential expression results. The results output a `gene column` along with the `logCPM`, `logFC` and `p-value` column. I have a question regarding conversion of log2 scale data to linear data for certain analysis. How can this be accomplished. Is there any R package to do this?

For instance, can I use the base function in R like below:

Linear FC value = 2^(logFC)

Linear CPM value = 2^(logCPM)

Thank you,

Toufiq

Answer: All results output from edgeR are on a log-2 scale, so yes you can unlog the logFC and logCPM values using the formula you give (as @rpolicastro already confirmed in his answer).

Let me say philosophically though that I don't agree that unlogged values can be described as "linear". Expression results are usually best analysed on a log scale and unlogged values do not behave in a linear fashion for most purposes.

Note it is especially important not to try to undo edgeR's logFC moderation, as mediated by the prior.count values, if you are working on an unlogged scale. If you really did want unmoderated fold-changes then you would simply set `prior.count=0` in the call to `exactTest` instead of doing your own ad hoc hacking. I can't imagine how that could be a good idea however. Getting infinite fold-changes from tiny counts (0 vs 1 for example) does not help anyone.
Question: Hello everyone,

I don't really know how to put this, but I'm recently trying to improve myself in coding and making tools. I was mostly doing a lot of bash things because it was the easiest way. But I know I should review my habits to be more portable.

This was my first year working on a cluster environment (SLURM), and I think I still don't get all the concepts behind it.
I was fine with just doing simple "sbatch", "srun" or "sarray" command, but I would like to understand more about how it works.

A concrete example :
I want to create a python script that call a tool (minimap2) installed on our cluster. This tool can only be used after a "module load" of the appropriate module (in that case bioinfo/minimap2-2.5).

So I was giving a shot with `subprocess.call` and `subprocess.run` (from which I understood they were sort of the same but not from the same python version ?) to perfom that module load.

But it kept failing, and I could'nt find anything about it on the net. I feel like I'm doing something abberant...

I'm not even sure if what I am trying to do is relevant (slurm module load from within python ?).
Can anyone correct me and maybe explain what I'm doing wrong ? Maybe if you have any website or book to recommand to understand better this topic. I feel very nooby, but I still want to improve.

Thanks for your help, sorry if not a proper Biostar question,

Roxane

Answer: The reason you are having trouble loading a module on the HPC with `subprocess.call`, etc., is because each subprocess is executed in a new, independent process. So even if you did a `subprocess.call` for something like "module load xyz", the loaded module will not persist into your following calls to use that program. If you really wanted to do that, you would have to do something like `subprocess.call('module load xyz; run_my_program')`, probably with the `shell=True` option enabled I believe.

Honestly your best bet is to use a workflow manager like Snakemake or [Nextflow](https://www.nextflow.io/). You should not couple your scripts directly to the HPC environment you are using, because they will then become unusable on other environments. For example on Nextflow, you can easily create a pipeline to run your script on your samples, and then configure Nextflow to submit the task as a batch job on [SLURM](https://www.nextflow.io/docs/latest/executor.html), and also configure it to first load the correct [modules](https://www.nextflow.io/docs/latest/process.html#module). Something like this would be the most portable, since it will allow you to update the execution configuration for a new system without having to modify any aspects of your pipeline's tasks.
Question: I have few motifs whose location on a given DNA sequence has been identified using MATCH (from TRANSFAC). MATCH provides the location of these matrices in a numerical way as motif 1 236 (+) (for an example - motif 1 is present at location 236). I want to get a graphical output by feeding the length of DNA sequence and the the position of matrices so that its easy for presentation and comparing two different sequences with almost similar motifs

Thanks in advance for any suggestions

Answer: Perhaps look at WebLogo, which assists with making sequence logos.
Question: I am dealing with a bunch of inversions and running some analysis on each of them .One of them is located on chrommose one, the starting position is 6199548 and ending position is 9699260. So the length of this position is 3499712 bp. I want to generate non-overlapping windows based on the length of this focal positions (3499712) and assign the focal position to one of these windows that overlap with the starting and ending position of my focal position. I am achieving this by using `GenomicRanges::tileGenome` function in `r`:

seqlength <-c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)
bin_length<-c(3499712)
bins<-tileGenome(seqlength,tilewidth=bin_length,cut.last.tile.in.chrom
= T)
#bins_u<-unlist(bins)
write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names
= FALSE, row.names = TRUE) }

This generates non-overlapping 3499712 bps windows across the genome.

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712
"2" "chr1" 3499713 6999424 3499712
"3" "chr1" 6999425 10499136 3499712
"4" "chr1" 10499137 13998848 3499712
"5" "chr1" 13998849 17498560 3499712

However I have a major problem here that I have not been able to sort it out. I want one of these windows exactly overlaps with the starting and ending position of my focal position (6199548:9699260), only the focal position not any other position outside this region. This means 2 windows (one before and one after the focal position) will have a different length which my analysis is fine with it. I just do not know how to do it?

This is my desired output:

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712 "*"
"2" "chr1" 3499713 6199547 2699834 "*"
"3" "chr1" 6199548	9699260 3499712 "*"
"4" "chr1" 9699261 13998848 4299587 "*"
"5" "chr1" 13998849 17498560 3499712 "*"

As you see the the second window is now smaller, the third window accommodates the focal position and the forth window is a little larges than bin size.
Does anyone have any idea how can I do this task in R? Help me please, I am really stuck!

Answer: I'm not sure how you did your calculation but for me it looks like your region is **3499713** bp in length.

site <- GRanges("chr1:6199548-9699260")

seqlength <- c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)

bin_length <- width(site)

## check length of region
print(bin_length)
##3499713

bins <- tileGenome(seqlength, tilewidth=bin_length, cut.last.tile.in.chrom=TRUE)

## extract overlapping bins
flank.bins <- subsetByOverlaps(bins,site)

## find new coordinates excluding region of interest
new.flank.bins <- setdiff(reduce(flank.bins),site)

## combine new bins with your region of interest
new.bins <- c(new.flank.bins,site)

## remove old bins
bins <- subsetByOverlaps(bins,new.bins,invert=T)

## add new bins
bins <- sort(c(bins,new.bins))

## check that it makes sense
##table(width(bins[seqnames(bins) == 'chr1']))

write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names= FALSE, row.names = TRUE)
Question: Hi, apologies for this basic question as I am new to the field. I have been checking my NGS data using FastQC and the checks fail on the Kmer content section. There appears to be a sequence TAGATCGGAA at position 90-100 bp in the reads which is enriched around 12 fold (Obs/Exp Max). However nothing shows up in the 'Overrepresented sequences' or 'Adapter content' sections of the report (complete flat line). The sequencing was done as BGI and I do not know what primers were used. My question is should I be trying to remove this Kmer sequence?

If I use grep on the first million reads in my fastq file to look for the sequence I only find 36 which appears to be quite a low number?

$ gunzip -c data1.fq.gz | head -4000000 | grep TAGATCGGAA | wc -l

Thanks for the help

Answer: FastQC often displays failed checks - it does not mean your data is bad, they are just warnings.

The only thing you should always really check is the quality along the length of your reads, which might force you to do some trimming. In this case, I would say that your reads are ok as they are, and you can proceed to the alignment phase of your workflow.
Question: The plot is from this publication, see page 2 the bottom of the plot with ACTG letters in different colours and sizes.

Could anyone clarify in layman's terms what it represents?

More importantly, how do we plot it using R or any other tools?

Here is the plot, if you can't access the paper:

Found related post at SO: Plotting a “sequence logo” using ggplot2?

Answer: Here is [**motiflogo**][1], which can make a SNP-specific motif logo representation. See figure, please.

[1]: https://github.com/zhilongjia/motiflogo
Question: Dear all,

I have a fasta file contatining marker ID's:

>marker1
AGCTGGGTGTCGTTGATCTTCAGGTCCTTCTGGATGTACAGCGACGCTCCActcatcacg
ccgatgtcgtagccgaggaggatggaggtcatggaggcgag
>marker2
CTCGTCGGTTTCGACCACCATGTTGTCATGGGCATAGCATAGTGAACTCCACTCACAGCT
AAAAAGGAACGGTGCCTGCCTTCTTTTCAGTATCTAATTAC
>markern
CATCTAATTTGTTGGACGGGGATATGCTGGGAAAACGAAGAGCAAATCTGTGCCCACTGT
GCTGGTCGCTGCAGAATCTATTCCCCATTATGTTGCCTGTT

I would like to align these markers against the reference genome and find the position of variant for every single marker if there any:

#CHROM	POS	ID	REF	ALT
chr1	1101439	marker1	C	A
chr1 1104710	marker2	A	T
chr1	1104748	markern	C	T

I am trying to use both `bwa/bowtie` to align and `bcftools/gatk` to call variants, but I do not know how to include also the marker ID in the vcf file?
OR MAY BE THIS IS POSSIBLE BY PARSING THE BLAST OUTPUT?

Answer: if you want to add the markerID, it's better to align with Blast and parse the result to generate the custom VCF
Question: Hi,

Is there a way I can introduce mutation in bam file at specific position of a chromosome without running the alignment. For example I would like to change base from A>C at Chr1:12345?

Answer: I quickly wrote: http://lindenb.github.io/jvarkit/Biostar404363.html

```
$ samtools view src/test/resources/toy.bam | tail -1
x6	0	ref2	14	30	23M	*	0	0	TAATTAAGTCTACAGAGCAACTA	???????????????????????	RG:Z:gid1

$ cat jeter.vcf
ref2	14	A	C

$ java -jar dist/biostar404363.jar -p jeter.vcf src/test/resources/toy.bam | tail -1
x6	0	ref2	14	30	23M	*	0	0	CAATTAAGTCTACAGAGCAACTA	???????????????????????	PG:Z:0	RG:Z:gid1	NM:i:1
```
Question: Hi,

I have managed to merge two different files using a common column which has the same information using the below function:

Combine_Files <- merge(File1, File2, by="Symbols")

However, I would like to know if its possible to merge two different *.csv files i.e. one file (file 1) with the gene symbols only (one gene per row), and the other file (file 2) with the gene symbols and other annotation columns (many related genes per row separated by comma but one of the gene is common as (file 1). I would like to merge both the files based on the mapping of the gene symbols from file1 and extract other annotations columns from (file 2) in the combined file for further data analysis. I would like to know how this could be done.

For example: File_1

Gene_Symbols

1. GeneA1
2. GeneA2
3. GeneA3

For example: File_2

Gene_Symbols

1. GeneA1, GeneX1, GeneX2, GeneD
2. GeneA2, GeneL1, GeneP2, GeneNA
3. GeneB3, GeneA3, GeneLP1, GeneNA1

Other columns in File_2

Phenotype, GO Ontology, Pathways

Expected Output

Gene_Symbols, Phenotype, GO Ontology, Pathways

Thank you,

Toufiq

Answer: Use the `tidyr` package. The function `separate_rows` will allow you to take you file_2 and split each row so that each row has only one gene symbol, and the information in the other columns is copied. Thus

GeneA1,GeneX1,GeneX2 WO cell processes Pathway 1

becomes

GeneA1 WO cell processes Pathway 1
GeneX1 WO cell processes Pathway 1
GeneX2 WO cell processes Pathway 1

You can then do your merge as before.
Question: Hello guys,

It may seem like a basic question, but this is causing confusion.

What is the difference between fold change and Log fold change?

Regards,

Answer: I find pictures helpful when thinking about these concepts, so quickly made some visualizations for you. They expand on WouterDeCoster's and h.mon's comments. The figures were made in R using with the following code:

#Generate fake data
fake <- data.frame("gene" = as.numeric(c(1:50)), "pre" = as.numeric(seq(1,99,2)), "post"= as.numeric(seq(99,1,-2)))
#Calculate fold change
fake$fold <- fake$post/fake$pre
#Calculate log2 fold change
fake$logfold <- log2(fake$fold)
#Visualize
library(ggplot2)
ggplot(fake, aes(x = gene, y = fold)) + geom_line()
ggplot(fake, aes(x = gene, y = logfold)) + geom_line() + ylab("Log2(FC)")

Let's say you have 50 genes (for the sake of convenience, the genes are called "1", "2", "3", ... , "50") and you measure expression before and after some type of treatment. You then measure the fold change of the genes due to treatment: FC = expression(post) / expression(pre). As russhh indicated, the logFC is simply log2(FC).

In this dataset, half the genes were upregulated and the other half downregulated. Interpreting the untransformed fold change is tricky: here it looks like gene1 had a huge fold change, close to 100, but what about genes30-50? It's hard to tell what their value is, and also by how much they differ from each other.

![enter image description here][1]

Performing the log2 transformation scales the data and facilitates the interpretation, as illustrated below.

![enter image description here][2]

The differences in expression in genes30-50 are now much more obvious.

[1]: https://i.imgur.com/dD7RC7n.png
[2]: https://i.imgur.com/5GywhTT.png
Question: Hi Biostar,

Apologies for this basic question but I can not find a simple answer

I want to sequence several samples of mRNA paired end 2x150 and for each sample I want 100 million reads (100 million for forward and 100 million for reverse)

the supplier proposes to sequence 3 samples per lane of flow cell of the Hiseq 4000 technology.
So what I understood 3 samples x 200 million reads (100 million for forward and 100 million for reverse) = 600 million reads per lane

but when I looked for the performances Illumina Hiseq 4000 I found this technology can generate 312000000 Clusters per lane only. How can this technology generate 600 million reads per lane ??

Would it be possible to give me a simple answer Thank you in advance

Answer: What is sequenced is a so called "fragment". To simplify, they will sequence each sample with 100 mln fragments, and each fragment is sequenced from both ends, i.e. paired end sequencing. And for each fragment you will have two reads - forward and reverse. So what the provider tells you is correct. Anyway, always ask them for details and explanations if you have doubts to avoid potential misunderstandings and waste of resources.