Reference file for RNA Seq Analysis

[PurplePancake]

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!

[dpryan]

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.

[PurplePancake]

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]

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>BBDDBCEC@B>@>CADA9B@ECB@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 GGEDGCEBEEHGHEBIIEIHDIFHIIGHIHHAFHIEIIIIGIIIFIHHIHIHIIIIIHIIIIIGIIIIIIIIIIII@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 HHHHHHHHHHHHHHDHHHHDDGGDGHHGHFHFFFHGDGDFHHHHHGH@CGEBHFEGDHEBBBDEC>B7CBBBB@B64*9==9>6==@=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<DI@IGFIIBDHIGFGGBDEGEB@GE<EDGDBG@GGGDEGDEHFDGIGE@BDGGBDGHGFIIIGG@II 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>BG@BGB>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]

One other detail

And the species is coral (A. digitifera)

[dpryan]

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.

[PurplePancake]

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!