cannin/biostars_qa huggingface dataset:
Question: I am trying to call peaks in ATAC-seq data. Not surprisingly, MACS is a popular option. According to the [MACS documentation][1]:

> 1. To find enriched cutting sites such as some DNAse-Seq datasets. In this case, all 5' ends of sequenced reads should be extended in both
> direction to smooth the pileup signals. If the wanted smoothing window
> is 200bps, then use '--nomodel --shift -100 --extsize 200'.
>
> 2. For certain nucleosome-seq data, we need to pileup the centers of nucleosomes using a half-nucleosome size for wavelet analysis (e.g.
> NPS algorithm). Since the DNA wrapped on nucleosome is about 147bps,
> this option can be used: '--nomodel --shift 37 --extsize 73'.

Based on a brief literature review, people use both `--shift -100 --extsize 200` and `--shift 37 --extsize 73` for ATAC-seq. Is one option more appropriate? Are there maybe different sub-types of ATAC-seq where one is better than the other?

[1]: https://github.com/taoliu/MACS

Answer: I got another alternative from the MACS developer:

> If you followed original protocol for ATAC-Seq, you should get
> Paired-End reads. If so, I would suggest you just use "--format BAMPE"
> to let MACS2 pileup the whole fragments in general. But if you want to
> focus on looking for where the 'cutting sites' are, then '--nomodel
> --shift -100 --extsize 200' should work.

And then more recently:

> Why do so many papers use --shift -100 --extsize 200 for MACS2 rather
> than -f BAMPE if this is not recommended for paired-end data?
>
> --shift -100 --extsize 200 will amplify the 'cutting sites' enrichment from ATAC-seq data. So in the end, the 'peak' is where Tn5 transposase
> likes to attack. The fact is that, although many information such as
> the insertion length and the other mate alignment is ignored, such
> result is still usable. Especially when the short fragment population
> is extremely dominant, the final output won't be off much.

Source: https://github.com/taoliu/MACS/issues/145
Question: I've been searching for some time, but I cannot find any md5sums that come with FASTQ files from the SRA. I'm currently using sra-tools' fasq-dump to get them from a script. On their website, only analyses seem to have md5sums.

Answer: The SRA archive format ("vdb") contains an md5 checksum as well as a few other consistency checks (I think). The sra-toolkit has a utility, [vdb-validate][1] which will report any errors in the data, and perform an md5 checksum comparison.

[1]: http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=toolkit_doc&f=vdb-validate
Question: Hello Everyone,

I have a fastq file and I want to extract only those reads which have length greater than 25 bp. So i want to make another fastq file with read length > 25 bp. How can I do this. This is my top 100 lines of fastq file

```
@SRR1024131.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1911:1053 length=100
AGGGCAAGTATGAAGAAGTAGAATATT
+SRR1024131.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1911:1053 length=100
DDFHHFHHGGGHGGIFHIJIIDIIJJI
@SRR1024131.2 DBRHHJN1:259:D0PM7ACXX:1:1101:2522:1198 length=100
GGCTCAACTTTCGATGGT
+SRR1024131.2 DBRHHJN1:259:D0PM7ACXX:1:1101:2522:1198 length=100
FFFGHHHHJJJJJGFIJF
@SRR1024131.3 DBRHHJN1:259:D0PM7ACXX:1:1101:3117:1165 length=100
ACATTTTTGAGTGCTTACTACAGT
+SRR1024131.3 DBRHHJN1:259:D0PM7ACXX:1:1101:3117:1165 length=100
FFFHHHHHHHIHEHHFGHFHHGII
@SRR1024131.4 DBRHHJN1:259:D0PM7ACXX:1:1101:3474:1075 length=100
TAGTACTTAGCAAAGAGTGA
+SRR1024131.4 DBRHHJN1:259:D0PM7ACXX:1:1101:3474:1075 length=100
DDDFHDFHIAGHIGHG@33A
@SRR1024131.5 DBRHHJN1:259:D0PM7ACXX:1:1101:3952:1099 length=100
TGAGAACTGAATTCCATAGGCTGT
+SRR1024131.5 DBRHHJN1:259:D0PM7ACXX:1:1101:3952:1099 length=100
EFFHGHHHHJIJJJJJIBFHEHIG
@SRR1024131.9 DBRHHJN1:259:D0PM7ACXX:1:1101:5277:1092 length=100
GCGGCGGCGTTATTCCCATGACCCGCCGG
+SRR1024131.9 DBRHHJN1:259:D0PM7ACXX:1:1101:5277:1092 length=100
FDDDHHDHI@B>=B>?@BD>ACCCBC@BB
@SRR1024131.11 DBRHHJN1:259:D0PM7ACXX:1:1101:6019:1101 length=100
AGTAGATTTGTATGGATTT
+SRR1024131.11 DBRHHJN1:259:D0PM7ACXX:1:1101:6019:1101 length=100
DDDHHFFHIGHAGHEFIII
@SRR1024131.14 DBRHHJN1:259:D0PM7ACXX:1:1101:8423:1248 length=100
AGTCGGTGATGGGAGTCTCT
+SRR1024131.14 DBRHHJN1:259:D0PM7ACXX:1:1101:8423:1248 length=100
FFFHHHFHIJIIJIJBHIJJ
@SRR1024131.15 DBRHHJN1:259:D0PM7ACXX:1:1101:9484:1233 length=100
TGCTGGGTCACACCTGAAGCT
+SRR1024131.15 DBRHHJN1:259:D0PM7ACXX:1:1101:9484:1233 length=100
FFFHHGHFHIJHHHJJHHIJJ
@SRR1024131.16 DBRHHJN1:259:D0PM7ACXX:1:1101:9807:1100 length=100
ACTATTCCAGCGAGAGTTAACATAAATTCCAAT
+SRR1024131.16 DBRHHJN1:259:D0PM7ACXX:1:1101:9807:1100 length=100
FFFHHHHHJJIJJJJIHHGHIJJGJJJJIIJJI
@SRR1024131.17 DBRHHJN1:259:D0PM7ACXX:1:1101:10857:1034 length=100
TAATCATTTTAATTGTACAGTTCAGTAATGT
+SRR1024131.17 DBRHHJN1:259:D0PM7ACXX:1:1101:10857:1034 length=100
B?3CDFBFFFFFIIF:EFHAHIC?FE+ABHH
@SRR1024131.19 DBRHHJN1:259:D0PM7ACXX:1:1101:13257:1082 length=100
ATGTGTTTGTAGGTTGTTTGTTGTCTTTA
+SRR1024131.19 DBRHHJN1:259:D0PM7ACXX:1:1101:13257:1082 length=100
DFFHHHHHJFHHIHHJFGHIJJIFIIIIG
@SRR1024131.20 DBRHHJN1:259:D0PM7ACXX:1:1101:14103:1161 length=100
TGAGGTAGTAGGTTGTATAGTT
+SRR1024131.20 DBRHHJN1:259:D0PM7ACXX:1:1101:14103:1161 length=100
FFEHFCFHFGHGEFHC@

Answer: I believe you could also use [seqtk][1] to do this:

seqtk seq -L 25 yourseqs.fastq.gz > cleanseqs.fastq.gz

[1]: https://github.com/lh3/seqtk/
Question: I am tying to retrieve mouse mm10 gene information using biomart library in R, but I don't know how to do that

(The information that I need are `mm10.knownGene.name`, `mm10.knownGene.chrom`, `mm10.knownGene.strand`, `mm10.knownGene.txStart`, `mm10.knownGene.txEnd` and `mm10.kgXref.geneSymbol`)

```r
source("http://bioconductor.org/biocLite.R")
biocLite("biomaRt")

library(biomaRt)

mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
listFilters(mouse)

getBM(attributes=c("ensembl_gene_id", "mgi_symbol"), filters= "mgi_symbol", mart=mouse)
```

Answer: If you just want mm10 symbol, chr, strand, transcript start & end, you could do this:

```r
res <- getBM(attributes = c("ensembl_gene_id", "mgi_symbol","chromosome_name",'strand','transcript_start','transcript_end'), mart = mouse)
```

If you have list of genes:

```r
#genesym is a character vector of gene symbols
res <- getBM(attributes = c("ensembl_gene_id", "mgi_symbol","chromosome_name",'strand','transcript_start','transcript_end'), filters = genesym, mart = mouse)
```
Question: I have in excess of 400 BED files of transcription factor binding site coordinates that I want to compare with one master BED file of binding regions. The aim is to identify which of the ~400 TFBS overlap with the master file.

Normally I would go straight to bedtools, but given the number of TFBS files I wonder if there is a "better" method? [bedtools intersect][1] would appear to do the trick, but for 400 files....?

The type of output I am looking to get is:

```
master_pos1 TFBS8, TFBS16, TFBS200
master_pos2 TFBS1, TFBS333
..
```

Thank you

[1]: http://bedtools.readthedocs.org/en/latest/content/tools/intersect.html

Answer: If your BED files are sorted, you could use [*bedops*][1] to union all the TFBS to standard output, and pipe that result to [*bedmap*][2] to do the mapping of TFs to master regions.

(This technique assumes that the TFBS files are minimally BED4. That is, the fourth column in each TFBS file contains the ID of the TF. If that is not the case, describe your format in more detail and I'll suggest a quick one-liner with *awk* to fix up files into the correct form.)

Here's the one-liner that unions and maps:

$ bedops --everything tfbs*.bed | bedmap --echo --echo-map-id-uniq --delim '\t' master.bed - > answer.bed

Piping to standard output avoids the unnecessary step of making an intermediate file somewhere on the hard drive, which is otherwise very expensive in time. So this should be very fast.

Assuming that the ID fields in each of `tfbs001.bed` through `tfbs400.bed` contain the desired TF names or other identifiers of choice, the file `answer.bed` contains the results as you expect, except that it uses a semi-colon as an ID delimiter, instead of a comma. You could add `--multidelim ','` to the *bedmap* statement, if that is a requirement.

If your BED files are not sorted, you could first prepare them with BEDOPS [*sort-bed*][3], which is faster at sorting BED files than GNU *sort*.

```
$ for tfbs_fn in `ls tfbs*.bed`; do sort-bed $tfbs_fn > sorted.$tfbs_fn; done
$ sort-bed master.bed > sorted.master.bed
```

Then use the sorted files in downstream BED ops. You only need to sort once.

[1]: http://bedops.readthedocs.org/en/latest/content/reference/set-operations/bedops.html
[2]: http://bedops.readthedocs.org/en/latest/content/reference/statistics/bedmap.html
[3]: http://bedops.readthedocs.org/en/latest/content/reference/file-management/sorting/sort-bed.html
Question: Dear all,

I am trying to generate exactly the same results in Limma for my result table as I did in DESEQ2 (Differential Expression Analysis).

With the toptable function I am not getting lfcSE and basemean.

Can anyone help me how to extract this information via Limma?

I do get logFC, AveExpr, t, P-Value, adj.Pvalue and B in Limma.

Thank you,
Bine

Answer: It would be superfluous to include a SE column in the limma output because it is not needed and because it is immediately computable from the other columns.
Question: hi,

I am trying to read my CEL files from [Mouse430_2] Affymetrix Mouse Genome 430 2.0 Array but I get this error

library(affy)

> Data<-ReadAffy()

Error in affyio::read_abatch(filenames, rm.mask, rm.outliers, rm.extra, :
Cel file C:/Users/Lenovo/Desktop/GSE50833_RAW/GSE10000_RAW/GSM44660.CEL/GSM252007.CEL does not seem to have the correct dimensions

> celfiles<- list.files("GSE10000/CEL", full = TRUE)

> rawData<- read.celfiles(celfiles)

All the CEL files must be of the same type.
Error: checkChipTypes(filenames, verbose, "affymetrix", TRUE) is not TRUE

how to read my CEL files for normalization?

Answer: The code you show is intriguing because you have in the error message a reference to the GSE10000 dataset, which is indeed Mouse430_2 arrays, but also to GSE50833, which are Agilent-028005 SurePrint G3 arrays. At any rate, the error suggests that you are trying to read different arrays with ReadAffy() and that fails because different arrays have different dimensions. First thing I would try myself is to make sure that all the files are of the same platform/array.

**EDIT**

I could replicate the problem and confirm my guess using the following experiment (there must be a better way to do this than reading the whole set of files one by one):

f <- list.files(pattern = "CEL.gz")
celf <- lapply(f, function(x) ReadAffy(filenames = x))
table(sapply(celf, annotation))
mouse4302 mouse430a2
18 17

The solution is to read them separately.

**EDIT 2**

OK, this is the most effective (fast) way to check the chip type of a bunch of cel files:

library(affyio)
f <- list.files(pattern = "CEL.gz")
table(sapply(f, function(x) read.celfile.header(x)$cdfName))
Mouse430_2 Mouse430A_2
18 17

**EDIT 3**

And this is how you can use the information above to read the files in different batches:

ff <- split(f, sapply(f, function(x) read.celfile.header(x)$cdfName))
ff
$Mouse430_2
[1] "GSM250879.CEL.gz" "GSM250880.CEL.gz" "GSM250881.CEL.gz" "GSM250882.CEL.gz" "GSM250919.CEL.gz" "GSM250920.CEL.gz"
[7] "GSM250922.CEL.gz" "GSM250923.CEL.gz" "GSM250925.CEL.gz" "GSM250927.CEL.gz" "GSM250928.CEL.gz" "GSM250943.CEL.gz"
[13] "GSM44658.CEL.gz" "GSM44659.CEL.gz" "GSM44660.CEL.gz" "GSM44661.CEL.gz" "GSM44662.CEL.gz" "GSM44663.CEL.gz"

$Mouse430A_2
[1] "GSM252007.CEL.gz" "GSM252008.CEL.gz" "GSM252009.CEL.gz" "GSM252010.CEL.gz" "GSM252011.CEL.gz" "GSM252014.CEL.gz"
[7] "GSM252015.CEL.gz" "GSM252016.CEL.gz" "GSM252017.CEL.gz" "GSM252018.CEL.gz" "GSM252021.CEL.gz" "GSM252022.CEL.gz"
[13] "GSM252033.CEL.gz" "GSM252040.CEL.gz" "GSM252051.CEL.gz" "GSM252052.CEL.gz" "GSM252053.CEL.gz"

library(affy)
abatch1 <- ReadAffy(filenames = ff$Mouse430_2)
abatch2 <- ReadAffy(filenames = ff$Mouse430A_2)

And so on.
Question: Dear community,

I have a huge paired-end HiC dataset (BAM format) which I want to format like this way:

```
HWI-D00283:117:C5KKJANXX:2:1101:1139:77789 chr6 153338506 153338556 37 + chr6 153338031 153338081 37 -
HWI-D00283:117:C5KKJANXX:2:1101:1139:77856 chr6 149915169 149915219 37 - chr6 149914908 149914958 37 +
HWI-D00283:117:C5KKJANXX:2:1101:1139:79414 chr4 184474969 184475019 37 - chr4 184474811 184474861 37 +
HWI-D00283:117:C5KKJANXX:2:1101:1139:81280 chr6 153641723 153641773 37 - chr6 153641551 153641601 37 +
HWI-D00283:117:C5KKJANXX:2:1101:1139:81917 chr8 87070282 87070332 37 - chr8 87069851 87069901 37 +
HWI-D00283:117:C5KKJANXX:2:1101:1139:82575 chr17 56970884 56970934 37 - chr6 151400450 151400500 37 -
HWI-D00283:117:C5KKJANXX:2:1101:1139:86642 chr6 150043041 150043091 37 - chr6 150042915 150042965 37 +
```

This is an example which I obtained by first converting the BAM format to BED and separating each mate into different files and then with a AWK command joined the mates.

This is the awk command I used:

awk 'NR==FNR {h[$4] = $1"\t"$2"\t"$3"\t"$5"\t"$6; next} {OFS="\t"; print $4,$1,$2,$3,$5,$6,h[$4]}' mate1 mate2

This command worked fine with a small dataset (1M, 10M reads), but when I tried with 200M reads file, it crashes because memory reasons I suppose. Is there a way to efficiently join paired-end reads as I showed in my example?

Thanks!

Answer: extract the mapped reads F, sort on name:

samtools view -f 64 -F 3976 your.bam | cut -f LC_ALL=C sort -k1,1 > F.txt

extract the mapped read R, sort on name:

samtools view -f 128 -F 3976 your.bam | LC_ALL=C sort -k1,1 > R.txt

join both files and format the output with awk:

LC_ALL=C join -t ' ' -1 1 -2 1 F.txt R.txt | awk -f your.script > result.txt
Question: I have a large number of short patterns (5-10 "letters" or amino acids in this context). I want to search for these in a large database of "texts" (proteins). Are there any existing packages which provide efficient algorithms for this? Ideally they would be

Handle 140k patterns and 280k texts

Able to handle ambiguous residues such as 'J' for [IL] (this can be faked by duplicating some patterns)

python is preferred (I'm already using BioPython, but it is missing this afaik)

I'm currently using a naive O(n^2) search, but it's taking too long. I think this problem is generally solved using keyword trees (Aho-Corasick algorithm), but I'd rather not implement that myself unless I have to.

Update

I should have specified that for output I need both the text and the pattern which matched it. This prevents some of the otherwise great implementations (eg fgrep, ahocorasick) from solving my problem.

Answer: I ended up using the esmre library, as suggested by dariober in the comments. It dropped my run time from hours (with a brute force search) to seconds.

Full code is on github, but the important parts are quite simple:

import esmre

vocab = esmre.Index()
for word in patterns:
vocab.enter(word,word)
for text in texts:
matches = vocab.query(text)
for match in matches:
print "Match for %s in %s" % (match, text)
Question: Like this:

```
HETATM 1519 C1 AAB1 B 100 -12.826 12.835 34.863 0.50 14.18 C
HETATM 1520 C1 BAB1 B 100 -4.549 17.796 20.909 0.50 14.18 C
________________^
```

There is *only one* C1 atom in AB1, but now we have the A and B version of this. I guess that because of the X-ray crystallography method, which has diffraction. But I am not sure about that.

Answer: When interpreting the x-ray diffraction and its correspondent electron density map, sometimes the position the a particular atoms is not crystal clear (pun intended). The atom can be too mobile, the density particularly bad at in that region, etc.. it's also tied to the occupancy factor (the second to last column, 0.50 in your example) which is 1.00 if the atom is observed 100% of the time at one place, less if split between different locations. In these cases, you have multiple entries per location of the same atom, differentiated by the character at position 17.
Question: Hi,
I have been trying to run RNA seq analysis on some paired end data. I have aligned on HISAT2, and run Stringtie, Stringtie Merge and then Stringtie again. To do the analysis I am using:
grch38_tran.tar.gz - https://ccb.jhu.edu/software/hisat2/index.shtml
Homo_sapiens.GRCh38.84.gtf - ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz

My issue is that despite running stringtie again after merge to remove some of the MSTRGs, I am getting a large number of them in my data set. More alarmingly the MSTRGs that do exist represent the highest counts in my sample.HISAT2-2.1.0.aligned.sorted.StringTie.1.3.3.gene_count_matrix.

Number of each: 24801 mstrg / 33970 ensg

Fraction of total: .42199 mstrg / .57800 ensg

Sum of each counts: 78615368 mstrg / 778402 ensg

Fraction of counts: .99019 mstrg / .00980 ensg

So while the MSTRG only makes up ~42% of the gene ids, it is 99% of what has been counted.
I have minimum coverage set to 5, and have -G set, as well as -e to restrict to the reference given.

Is there anyway to further optimize this? Have I missed out on an important step?

Answer: This has always been an issue as far as novel transcript discovery goes, you can see a lot of hits. Keep in mind that the vast majority of these are very slight changes to known transcripts and splice events, which are generally meaningless. When performing this kind of analysis I generally get rid of any MSTRG ID that falls within a known annotation, and then look for protein coding potential of transcripts identified, finally prioritising on abundance. I'll then go through a short list of these transcripts and visualise them in IGV to see if they're convincing.

A lot of this prioritisation I've been able to do with awk, and drastically reducing noise with [TACO][1], as a replacement for stringtie-merge. TACO also includes a utility to compare your merged GTF against a reference GTF, which is handy for subsetting.

[1]: http://tacorna.github.io/
Question: I like to know is there any way to merge the same coordinate in a table.

Table:

chr1	155944562	155945214	fantom_neuron	GSM1554667	9.84447
chr1	155944562	155945214	fantom_neuron	GSM1554672 7.43630
chr1	155944562	155945214	fantom_neuron	GSM1554678	32.77627
chr1	155945743	155946196	fantom_neuron	BAMPE	3.87072
chr1	155945743	155946196	fantom_neuron	GSM1554666	18.14939
chr1	155945746	155946202	fantom_neuron	GSM1554655	1.14939

Expected table:

chr1	155944562	155945214	fantom_neuron	GSM1554667,GSM1554672,GSM1554678	9.84447,7.43630,32.77627
chr1	155945743	155946196	fantom_neuron	BAMPE,GSM1554666	18.14939 ,3.87072
chr1	155945746	155946202	fantom_neuron	GSM1554655	1.14939

Answer: Using **R**, and *data.table* package:

library(data.table)

# use fread for your data
# mydata <- fread("myFile.bed")

# example data
mydata <- fread("
chr1 155944562 155945214 fantom_neuron GSM1554667 9.84447
chr1 155944562 155945214 fantom_neuron GSM1554672 7.43630
chr1 155944562 155945214 fantom_neuron GSM1554678 32.77627
chr1 155945743 155946196 fantom_neuron BAMPE 3.87072
chr1 155945743 155946196 fantom_neuron GSM1554666 18.14939
chr1 155945746 155946202 fantom_neuron GSM1554655 1.14939
")

# then group by paste
mydata[, lapply(.SD, toString), .SDcols = c(5:6), by = list(V1, V2, V3, V4) ]
# V1 V2 V3 V4 V5 V6
# 1: chr1 155944562 155945214 fantom_neuron GSM1554667, GSM1554672, GSM1554678 9.84447, 7.4363, 32.77627
# 2: chr1 155945743 155946196 fantom_neuron BAMPE, GSM1554666 3.87072, 18.14939
# 3: chr1 155945746 155946202 fantom_neuron GSM1554655 1.14939
Question: Hi I have a set of reads mapped to a reference sequence, what I what to do is to take the reads that map to one sequence and assemble them in a comtig as best as possible. Then I want to repeat this for thousands of sequence-read sets. I do not want to put this in a regular assember like mira, abyss, etc. because it is an overkill for what I want to do. I also do not want any software messing with my assignment and going its own way. Is there such a command line tool?

Thanks

Answer: TASR and SSAKE are traditional bets for this sort of thing, as well as PRICE. Haven't tried Tadpole, as genomax2 suggests, but I'm sure that's worth a look as well - especially as ones I've mentioned aren't going to break any speed records!
Question: Hi there
Can anyone explain to me how to use the ESTIMATE package in RNA-seq analysis?
I want to calculate immune scores and stromal scores by employing the ESTIMATE algorithm, then analyze the relationship of immune/stromal scores with subtype classification and cytogenetic risk by one-way analysis of variance, but I don't know how to do this!

I will be grateful for any help you can provide.

Answer: A vignette PDF comes installed with the package, and should be located at:

- R/x86_64-pc-linux-gnu-library/4.0/estimate/doc/ESTIMATE_Vignette.pdf

In this vignette, they use some data that comes bundled with the package (*R/x86_64-pc-linux-gnu-library/4.0/estimate/extdata/sample_input.txt*), which represents Affymetrix U133 microarray data that appears to be normalised and transformed by log [base 2].

So, if you have RNA-seq data, I would normalise the data in the usual way, and then transform via rlog or vst. Then, with ESTIMATE, use the rlog or vst expression levels.

Kevin
Question: I have several CPU nodes connected to the same network mounted directory. I want to parallelize the abyss-map step after abyss generates unitigs.

When I run abyss-map on my reads, it outputs a message saying it is generating a suffix array for the mapping. Does it store that on disk? Or does it just store it in memory for the mapping?

If I parallel run several abyss-map in the same directory, will the suffix array generated overwrite the suffix array generated by other abyss-map runs?

Answer: Hi Damian,

`abyss-map` builds the indexes in memory by default, so there should be no conflict about running multiple instances on the same file in parallel.

However, you can also use `abyss-index` to pre-build the indexes (.fai and .fm files). Any subsequent `abyss-map` runs will re-use those index files (similar to how it works with BWA), which allows you to run `abyss-map` more quickly and with a smaller memory requirement. (The initial indexing step requires an amount of RAM that is about 10 times the size of the input FASTA file.)
Question: Hello,

I am trying to analyze the public dataset https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126030
I've downloaded the fastq files onto my cluster, and would like to proceed with cellranger count.

I am in a test folder and the only file is: SRR8526547_1.fastq and refdata-cellranger-GRCh38-1.2.0

cellranger count --id=cellranger \
--transcriptome=/home/jl2/scratch60/refdata-cellranger-GRCh38-1.2.0/ \
--fastqs=.\
--sample=SRR8526547_1.fastq \

I keep getting the error of

Invalid path/prefix combination: /gpfs/ycga/scratch60/k/jl2/test, ['SRR8526547_1.fastq']
No input FASTQs were found for the requested parameters.

Can't seem to figure out what's wrong. Does it need fastq.gz instead of fastq?

Answer: `cellranger count` expects a certain nomenclature for the `fastq` files, please see the last section [here](https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input), "My FASTQs are not named like any of the above examples".

Basically this is how your file names should look like: `[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz`.

For the `Read Type`, you can take a look at your fastq files with `head` to see what is what. The link above explains different read types.