SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Downstream RNA-seq analysis without reference genome nareshvasani De novo discovery 26 05-21-2014 06:23 AM
Can I use a small subset of transcripts as a reference for RNA-seq analysis JackMetal RNA Sequencing 4 09-25-2013 06:47 AM

Reply
 
Thread Tools
Old 08-09-2014, 12:10 AM   #1
PurplePancake
Member
 
Location: New York City

Join Date: Aug 2014
Posts: 14
Default Reference file for RNA Seq Analysis

Hello,

I have .sam files that I would like to feed into htseq-count. However, I must first have a reference file in .gtf format, I believe.

I have access to an "assembled transcriptome" that is .fa.gz format of the same species from the .sam files. Could this possibly be used for the htseq-count as reference sequence input, and if so, should I do something to change its format as an acceptable input?

I am sorry for the noobie question!
PurplePancake is offline   Reply With Quote
Old 08-09-2014, 02:53 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

Did you just map to the assembled transcriptome? If so, you don't need to use htseq-count and can instead just use "samtools idxstats". You might need to filter by MAPQ or something to get rid of multimappers, but that'll be simpler (you'll just output this to a new file and index that).

If you didn't map to the assemble transcriptome but instead mapped to the genome, then we'll need some more details about your experiment (in particular, the species used). htseq-count needs an annotation file because it was designed to be used for experiments where reads are mapped to the genome.
dpryan is offline   Reply With Quote
Old 08-09-2014, 06:04 AM   #3
PurplePancake
Member
 
Location: New York City

Join Date: Aug 2014
Posts: 14
Default Thanks you!

Hello, I am not sure what was done to the file (I downloaded it). Is there some way I could determine that? The first 50 lines are in the format:

>adi_EST_assem_1
AGTAACACTGTGAAGTTGTAAAAATTCTGCATGAAAAAGCAGTGTTGTGACTCATAGTAT
GCAAAGAACTTTCGCGAGAATCTTTTCCACAAGTTGTTGTGACTAGATTTGAAAGAGAAG
AATGGGTATGACGTCAGTTTCAGGCGATCTACCTGATATTATGTGGAAATCAAGATGTTG
CTGGGAGTCCCAAAATGAATACAAACATAACAATAGCTTGTCATGCCTGAGTTATTCCAA
AGCATATGCGAACGGCAAATGGAACTCGGTAAAATGAAGCTTATCTCTTGCATGTTATAC
CTTCTACATAGCGTCGTATTATTAGATGGATTTCAGTGGTCTCATCGACAACAGGCAAGC
CACGTTATCAGCGTTGAAGCTGGTTCAGATGTTTCTCTTCCATGCGAATACGAACTCACG
CCACAAGAACAGCAGGAAGCTGATGTCTTTCATCTATTAACATGGACACGAGAAGAGCCT
CATTTCAGCGATCGGTGGACCGGATTGGCGATAAATTCGACTCTAAGTGAAAGCAAAGTT
ATTTATGATGATCCACAGCGAATTTTCATCGTTGATGGAACACTTACTGTCAAGAACATA
ACCGTTAAAGATAGGACTCGGTATCAATGCGCATTTCAAAGCAGTTTCTTCACGACACCC
AGCATAATTAACTTGGACGTCAAATATCCTCCAGAGATTATGCTCAAACCATTGTCCCAA
GATATCATTGAAGGGAATGATGTCAGGCACTGTGCAATGCATCTGGAAATCCACAGCCTA
AGATAACTTGGACGAGGCAAGGAAACAGTGATGTCCTTTCTTCATCTGAGACCCTTTTGT
TGCGTAATTTAGTCAGCGAAGATGATGGATCAGTGTATACGTGCAGGGTCGAAAATTACC
TTGGATTGAAGCAAGCCTCTGTTACAATAACCATGCAATATAAGCCAAAGGACACACGAG
TTACTATTTCATCTGGTGAGGTGAAGCTAGGAGATCGTGTGGACATCTCATGCTGGGCAC
GTGCAAACCCTTCACCAGAGTACAAATTTTATTATAACAACAAACTGATCAGGTGGTCAA
GTAAAGGACTGCTTTCATTAACGAGTGTTAATACTGAAGATCAAGGCACGTATCAATGTG
TTCCAGTGAATTACCTGGGGGATGGATCTGGGGCATCTGTTACTGTAACTTTGTCCACTG
GTGACAAGGCCTTCCCAGTATGGGCGTATGCTGCAATAGCAAGTGGTGTTCTCTTCCTTC
TTTTAATTGCTGCGGGAATTGCATTGTGCAGGGCAAGAAAAAAATGGAAGGCAGCAGATG
ACGGGACTAAGACAGCCCAACCGAAGGAGAGCACAATTAGCAGGCCTAGAGATGCTAAAG
AGGCACATGTAGGTCCAGGAGATTTCTATGGAGACCACAATGGCAGTATTGAGATTCCTC
ATCTTAATGGAGCTCTGAGTTTCAGTGACATCAGAAGTAGTCAAGAGTTCCCAGGTTTAA
AGGGTGCTTTT
>adi_EST_assem_2
AGTTCAGAACGGATTGCGTGGTTTTGGCAAAAGATAGAACAGCTGAGTTGTGAAAGGAAA
TTTGGAGTAAAGTAATTGGAATGGGTAGAAAGAACAGCAGATGTAGCCTAGCTGGGGTCG
TCAGCCTGAGCGTACTCTTCTGCATCGTAAGCTTTTTCATGATCGTTGACGCATCACATC
AAAGTAAACGAACACAGAGGCGAGAGAGTCTCGACAATAACAGAGTGAACTTGAGCGCCG
GAAATTTCCTCGTAGAATGGAACGTTGACAAAAACTCTGGAAGCATTGAATTCATATTGA
AGGTTACAGTTGATGAAAAAGGCTGGATTTTCTTTGGATTCATGCCCATCTCAATTAACA
GTTCCACGAAGCCAAGCGGCATTGAAGATTCTACAGGAGACTTCGTGGTAACATGGTTAT
TACCTGGCTCTAACGACAGAGAAACTTTGGATTTGAATACAATAAAAAAGAAAGGGGAGC
TAGGGATAGACACCGACAACAATTATGGGATTACAGCAGGGGATTCTAACAAGAATTCAA
GTGTACAGGTGCTTCACATTACCAGGAAACTTGACACAGGCGACCCACAAGATGTGCCCA
TCACAAACCAACCAATGTGGATGTTTGCTGCATGGGGAACATCTGAAGAATTCACAAAAA
ACTTCACCACAATTAAGGACAACATAACCAACTCCCTAGTGGTGGTAAAAGTGGTTTTTT

I'm guessing that does not help, though, right?

Where I downloaded the file, it said "Assembled Transcriptome v1.0, (compressed with gzip), total 36,780 contigs for 29,364,984 bp"

I am doing this just for practice to get a feel for RNASeq. But I have limited time, so I am not too concerned about quality at all (NOT to be published etc)!!
PurplePancake is offline   Reply With Quote
Old 08-09-2014, 08:42 AM   #4
PurplePancake
Member
 
Location: New York City

Join Date: Aug 2014
Posts: 14
Default Clarification

dpryan:

Thanks again.

I realize I did not show what the .sam files look like. The head is as such:

@HD VN:1.0 SO:coordinate
@SQ SN:scaf1 LN:42540
@SQ SN:scaf10 LN:28213
@SQ SN:scaf100 LN:68029
@SQ SN:scaf1000 LN:50671
@SQ SN:scaf10000 LN:2123

and tail is as such:

WUSI-EAS1631:100:FC63B34AAXX:1:99:14443:9257 0 scaf9999 4971 50 21M406N79M * 0 0 CAAAACGCTGACGCCATCCACCTCAACCAAATGGTCTCCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGC HHHHHHFHDHHHHHHHEHHHHHBHGHHHHHHHFHEDFFGFHBHFHCFED>[email protected]>@>[email protected]@BEAEBA<@@@@############### XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1
HWUSI-EAS1631:100:FC63B34AAXX:8:13:5760:14033 16 scaf9999 5406 50 87M508N13M * 0 0 AATGGTCTCCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTG GGEDGCEBEEHGH[email protected]IIIIIIIIIHIIIIIIIIIIIII XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1
HWUSI-EAS1631:100:FC63B34AAXX:2:68:12267:10621 0 scaf9999 5414 50 79M508N21M * 0 0 CCTGGCTCGACGACTCCATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCCTCTAATGAATTTTGAGTTCCTT [email protected]C>[email protected]*9==9>[email protected]=BBB==<===4;= XA:i:1 MD:Z:77A22 NM:i:1 XS:A:- NH:i:1
HWUSI-EAS1631:100:FC63B34AAXX:2:67:17937:17163 16 scaf9999 5424 50 69M508N31M * 0 0 CGACTCCATTCTGCCCCACCGGCCCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTGAGTTCCTTCCAAACTAAT ####################@?A><EC3GD<[email protected]@GE<[email protected]@[email protected] XA:i:2 MD:Z:13A8A77 NM:i:2 XS:A:- NH:i:1
HWUSI-EAS1631:100:FC63B34AAXX:1:99:14443:9257 16 scaf9999 5431 50 62M508N38M * 0 0 ATTCTGACCCACCGGACCCTCAGGACTCACGGCATGCACTCTGTGCCATGGCAAGGAACCATCTAATGAATTTTGAGTTCCTTCCAAACTAATTCCTAAG >B>?8G5DBFC<BFDFBIDE>[email protected]>BGGEDGCEE?C?EAGDHEDHIIIIHIGBGGGGDGG8GGGBGHDIIHGIFFDGGBGGBGGGGGGGGBGEGEGE XA:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1

Could that indicate whether the .sam files have been mapped? If so, then are you saying I can run "samtools idxstats" on the .sam files?

I did try to run samtools idxstats on them using the command:

../../samtools-0.1.7_i386-darwin/samtools idxstats adu.sam

But get the error:

[main] unrecognized command 'idxstats'

However, samtools did seem to be working when I earlier converted between files using samtools views, as such:
../samtools-0.1.7_i386-darwin/samtools view -h -o adu.sam adu.bam

I'm not sure what might be causing that error (yes, I am very new to this)...sorry (

Thanks so much again, for any help......!!!
PurplePancake is offline   Reply With Quote
Old 08-09-2014, 08:45 AM   #5
PurplePancake
Member
 
Location: New York City

Join Date: Aug 2014
Posts: 14
Default One other detail

And the species is coral (A. digitifera)
PurplePancake is offline   Reply With Quote
Old 08-09-2014, 09:57 AM   #6
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,479
Default

You can probably use the gff file. The genome sequence associated with it is here, in case it's different from the one you aligned against.
dpryan is offline   Reply With Quote
Old 08-09-2014, 01:53 PM   #7
PurplePancake
Member
 
Location: New York City

Join Date: Aug 2014
Posts: 14
Default Thank you!

dpryan:

Thank you so much!

I was wondering, where did you find this .gff file?

I followed steps to successfully install htseq-count. I believe it is actually working (rare for me), as I can type into the terminal htseq-count, and it auto-fills as a command.

However, I tried a simple command like:

htseq-count adu.sam aug_repeatmask_pasa_input.gff3 > ADUcount.out

and got the error:

Feature aug_v2a.20244.t1.exon1 does not contain a 'gene_id' attribute
[Exception type: ValueError, raised in count.py:53]

I am assuming the .gff3 file is missing a "gene_id" field, but am not sure how to fix that?

Sorry for another question, and thanks for your help with this!
PurplePancake is offline   Reply With Quote
Reply

Tags
htseq-count, rnaseq, transcriptome

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:42 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO