0
$\begingroup$

I did manage to bild index file (from NCBI transcriptome) and perform alingment in Rsubread in R with my fq files. I did get BAM files as a result and no error. However i am having trouble with next step as featurecount function does not produce any result( 0%). I have tried builtin annotation (hg38)

i am wondering if BAM files are ok..

this is the code i have used

library(Rsubread)
index_path <- "D:/FASTQ_Raw/hum_transcr"
all_files <- list.files(path = "D:/FASTQ_Raw", pattern = ".fq$", full.names = TRUE)
common_prefixes <- unique(gsub("_[12]\\.fq$", "", basename(all_files)))

output_directory <- "aligned_reads/"

dir.create(output_directory, showWarnings = FALSE)

for (prefix in common_prefixes) {
 
  read1_file <- paste0("D:/", prefix, "_1.fq")
  read2_file <- paste0("D:/", prefix, "_2.fq")
  
  output_file <- paste0(output_directory, prefix, "_aligned.bam")
  

  align(index = index_path,
        readfile1 = read1_file,
        readfile2 = read2_file,
        output_file = output_file,
        nthreads = 4) 
}

so the results/stats of individual BAMs appear ok still when i try featurecount function on BAM files i get zero....i was also using built in annotation index hg38 for attempt of counting but it does not work. I had issues with getting gtf file from NCBI. What do you recommend ?

$\endgroup$
4
  • $\begingroup$ How did you run featureCounts() and what was the output? Are you sure your reference and gtf file are compatible? $\endgroup$
    – Cloudberry
    Commented Jan 9 at 20:08
  • $\begingroup$ Well i have created Index file in first step from .fa transcriptome file from NCBI ( for gh38.14). Howeverr there i had problem of opening gtf file. I did download one from ensembl but i am not sure if they are 100 % compatible ( they should be i guess). Because of this i have tried Builtinannotation function from Rsubread but this yielded 0% alignment. I was hoping to do everything in R because i am beginer with seq analysis and with R I have at least some experience with different packages. $\endgroup$
    – MKE1508
    Commented Jan 9 at 20:38
  • $\begingroup$ NCBI and Ensembl use a different format for chromosome names. That's probably why you didn't get any results. The fasta and gtf must be from the same dataset. $\endgroup$
    – Cloudberry
    Commented Jan 9 at 20:46
  • $\begingroup$ Many Thanks!! I will retry with NCBI gtf file. Alignment was quite reasonable (over 90%) so it should work... $\endgroup$
    – MKE1508
    Commented Jan 9 at 20:53

1 Answer 1

-1
$\begingroup$

I don't think FeatureCounts is meant to work on an alignment to the transcriptome. It's not at all clear to me that the subread aligner, or any aligner, is meant to be used on a transcriptome.

If you have the transcriptome, I'd just use pseudoaligners like Kallisto or Salmon.

$\endgroup$
3
  • $\begingroup$ well i think it does, at least according to the user manual. And feature counting is mandatory step after alignment to the index. Not that i am an expert .IS Kallisto Windows compatible? I am not familiar with Linux... $\endgroup$
    – MKE1508
    Commented Jan 9 at 20:43
  • $\begingroup$ Aligners can definitely be used on a transcriptome. That's done all the time on non-model organisms for which full genome is not available. I don't think the Rsubread aligner is different in that sense. $\endgroup$
    – Cloudberry
    Commented Jan 9 at 20:55
  • $\begingroup$ But they aren't going to handle ambiguous alignments properly, like a psudoaligner would. $\endgroup$
    – swbarnes2
    Commented Jan 10 at 16:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.