參考這些
主要根據(jù)這篇操作！！！
??????我覺得學(xué)RNA-seq必看的文獻(xiàn)?。??????
這篇寫的特別好：RNAseq-workflow
RNA-seq原始數(shù)據(jù)的下載方法
 fastq格式大全
 本科生搞定RNA-seq上游數(shù)據(jù)分析
 RNA-seq一般流程
 轉(zhuǎn)錄組入門：序列比對
 RNA分析簡潔版
 序列比對

我覺得這 new_workflow/
│ └── annotation/ │
│ └── genome/ │
│ └── input/ │
│ └── output/ │ ├── 1_initial_qc/ │ ├── 2_trimmed_output/ │ ├── 3_rRNA/ │ ├── aligned/ │ ├── filtered/ │ ├── logs/ │ ├── 4_aligned_sequences/ │ ├── aligned_bam/ │ ├── aligned_logs/ │ ├── 5_final_counts/ │ ├── 6_multiQC/ │
│ └── sortmerna_db/ │ ├── index/ │ ├── rRNA_databases/ │
│ └── star_index/ 個作者這樣的分類特別清晰，值得學(xué)
<- Genome annotation file (.GTF/.GFF)
<- Host genome file (.FASTA)
<- Location of input RNAseq data
<- Data generated during processing steps
<- Main alignment files for each sample
<- Log from running STAR alignment step
<- STAR alignment counts output (for comparison with featureCounts)
<- Sequences that aligned to rRNA databases (rRNA contaminated)
<- Sequences with rRNA sequences removed (rRNA-free)
<- logs from running SortMeRNA
<- Main alignment files for each sample
<- Alignment files generated from STAR (.BAM)
<- Log from running STAR alignment step
<- Summarized gene counts across all samples
<- Overall report of logs for each step
<- Folder to store the rRNA databases for SortMeRNA
<- indexed versions of the rRNA sequences for faster alignment
<- rRNA sequences from bacteria, archea and eukaryotes
<- Folder to store the indexed genome files from STAR

??前半部分是我的血淚史！正文在中間開始！??

這里都是坑，別學(xué)?。。。≌_的在下一條分界線那里↓

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

其實序列就是文本文件（圖片來源于網(wǎng)絡(luò)）

1.2 轉(zhuǎn)為`fastq`格式

fastq-dump *.sra

生成一個fastq文件

到這里才完全發(fā)現(xiàn)問題?。。。?！??????

我的數(shù)據(jù)其實是有問題的！！！
我以為是單端測序（SE），但是在NCBI上查詢了一下原來是雙端測序（PE）?。。?br>

我選的是SRR3589959

所以就是重頭來過?。。?！??

?所以這里才是正文的開始！！

1 原始數(shù)據(jù)下載

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

SRR_Acc_List.txt里面是要下載的SRR號

注意！?。?！光標(biāo)一定不能在有SRR號的那行，不然會跳過那行??！

cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} );done

下載好了

1.2 `sra`轉(zhuǎn)成`fastq` (這里要注意不是這么寫?。。?！

因為是雙端測序數(shù)據(jù)?。?！

fastq-dump *.sra  絕對不能這么寫?。。。?

應(yīng)該加這個參數(shù) `--split-files`

fastq-dump --split-files *.sra

??????這里我會把.sra后綴的文件都移到一個新文件夾sra里，然后cd到sra文件夾里轉(zhuǎn)fastq↓

cp ./SRR*/.sra ./sra
fastq-dump --split-files *.sra

如果需要壓縮，可以加`--gzip`，就會生成`.fastq.gz`文件，解壓縮用`gunzip`

--split-files的意義

這樣就會是兩端，兩個文件

單端測序SE用`-U`參數(shù)

而且不知道是不是我的錯覺，好像分開兩個文件的速度比合并在一個文件中速度快一些

1.3 `md5sum`檢測傳輸?shù)?code>fastq文件是否完整，一般會附帶一個MD5校驗文件

 md5sum *.fastq | tee md5sum.txt

2 質(zhì)控

2.1 `fastqc`生成質(zhì)控報告 fastqc.html和 fastqc.zip格式

  fastqc -o . *.fastq

-o表示輸出文件的位置

會顯示這樣

怎么看結(jié)果

參考這篇
fastqc質(zhì)控

multiqc將各個樣本的質(zhì)控報告整合為一個

multiqc *.zip

生成一個文件夾一個html

multiqc的結(jié)果是兩個放在一起顯示

2.2 `trim_galore`修剪數(shù)據(jù),用于去除低質(zhì)量和接頭數(shù)據(jù)

參數(shù)--fastqc：在數(shù)據(jù)過濾后再次質(zhì)檢

#需要新建一個clean文件夾
mkdir ./clean
cd ./clean
ls 你放fastq數(shù)據(jù)的路徑/sra/*_1.fastq >1
ls 你放fastq數(shù)據(jù)的路徑/sra/*_2.fastq >2
paste 1 2  > config

如果是單端就直接↓
ls fastq數(shù)據(jù)位置/*.fastq > config
#config數(shù)據(jù)里其實就是fastq文件的路徑

#linux里輸
cat config |while read id
do
#二選一??
#雙端
 trim_galore ${id} -q 25 --phred33 --length 36 --stringency 3 --paired -o ./
#單端
 trim_galore ${id} -q 25 --phred33 --length 36 --stringency 3 -o ./
done 

#或者也可以寫成.sh，然后
bash qc.sh config

可以自動識別接頭

3 比對到參考基因組

我們有時糾結(jié)是否真的需要比對：
如果你只需要知道已知基因的表達(dá)情況，那么可以選擇alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）

打分機(jī)制：reads比對到基因組上的一個位置，我們稱之為一個alignment。軟件會對所有的alignments進(jìn)行打分和判斷，能夠符合過濾條件的alignment稱之為valid alignment，只有valid alignments, 才會輸出。
和BLAST類似，但在面對巨量的短序列數(shù)據(jù)時，類似BLAST這樣的軟件實在太慢了！因此，需要更加有效的數(shù)據(jù)結(jié)構(gòu)和相應(yīng)的算法來完成這個搜索定位的任務(wù)
每個alignment也有對應(yīng)的打分機(jī)制。hisat2 從以下幾個方面對alignment進(jìn)行打分
1 錯配堿基罰分
錯配堿基的罰分通過--mp參數(shù)指定，其值為逗號分隔的兩個數(shù)字，第一個數(shù)字為最大的罰分，第二個數(shù)字為最小的罰分
2 reads上的gap罰分
gap的罰分通過分成兩個部分，第一次出現(xiàn)gap的罰分和gap延伸的罰分，reads上的gap罰分通過--rdg參數(shù)指定，其值為逗號分隔的兩個數(shù)字，第一個數(shù)字為gap第一個位置的罰分，第二個數(shù)字為gap延伸的罰分。
3 reference上的gap罰分
reference上的gap罰分通過--rdg參數(shù)指定，其值為逗號分隔的兩個數(shù)字，第一個數(shù)字為gap第一個位置的罰分，第二個數(shù)字為gap延伸的罰分。
經(jīng)過一系列的罰分機(jī)制，每個alignment會有一個對應(yīng)的得分，然后會根據(jù)一個閾值，來判斷這個得分是否滿足valid alignment的要求。

hisat通過--score--min參數(shù)指定該閾值，指定方式是一個和reads程度相關(guān)的函數(shù)，默認(rèn)值為L,0,-0.2, 對應(yīng)函數(shù)為f(x) = 0 - 0.2 * x
根據(jù)reads長度，可以計算出得分的閾值，大于該閾值的alignment 被認(rèn)為是valid alignment , 才可能被輸出。L代表線性函數(shù)，此外，也支持其他類型的函數(shù)，比如常量，自然對數(shù)等，更多選擇請參考官方文檔。
一條reads可能會擁有多個valid alignments, 在輸出時，并不會輸出所有的alignments, 而是只輸出-k參數(shù)指定的N個alignments，-k參數(shù)的默認(rèn)值為5。
輸出結(jié)果以SAM格式保存，默認(rèn)輸出到屏幕上，可以通過-S參數(shù)指定輸出文件。

hisat2 -t -x genome路徑/genome -1 SRR路徑/SRR3589959_1.fastq.gz -2 SRR路徑/SRR3589959_2.fastq.gz -S 輸出路徑/SRR3589959.sam

-S表示輸出sam文件
只需要寫genome就好，它識別的是前綴，不是某一個文件！

hisat2 -t -x

但是報錯了，需要re-run on a computer with more memory

有人建議重新下載index，于是我在Hisat2官網(wǎng)找到了index
hg19(grch37) 和 hg38(grch38) 的區(qū)別

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grch38.tar.gz
tar -zxvf grch38.tar.gz

壓縮包有make_grch38_tran.sh文件，詳細(xì)記錄了創(chuàng)建索引的過程。

grch38

但是還是報錯，可能真的是內(nèi)存不足
未完待續(xù)……

再更新：
也有可能是index的問題？
可是我這個index是直接從hisat2官網(wǎng)下載的啊……

我決定先放棄Hisat2，試試bowtie2
建立矩陣

bowtie2-build 基因組路徑/hg19.fa index

但是建立不起來

時間太長了

試試`STAR`

參考這篇

STAR --runThreadN 6 --runMode genomeGenerate --genomeDir 你的路徑/hg19/ --genomeFastaFiles 你的路徑/.fa --sjdbGTFfile 你的路徑/.gtf --sjdbOverhang 100

--runThreadN ：設(shè)置線程數(shù)
--genomeDir：index輸出的路徑
--genomeFastaFiles ：參考基因組序列
--sjdbGTFfile ：參考基因組注釋文件
--sjdbOverhang ：這個是reads長度的最大值減1，默認(rèn)是100

應(yīng)該是內(nèi)存不足報錯了

這位朋友說加一個加內(nèi)存的參數(shù)
--limitGenomeGenerateRAM 80000000000（運行內(nèi)存加到80G）
--genomeChrBinNbits 11 這個的計算方式 : log2((基因組長度)120,000,000,000/(NumberOfReferences)10,253,694)=11
但是還是不行，又看到大家的討論，里面有這么一位網(wǎng)友↓

直到我看到了這個

所以我覺得在我的小小虛擬機(jī)里應(yīng)該是跑不動了，放棄 :-(

這篇作者在評論里建議用salmon
具體操作參考這篇

Salmon是Alignment-free transcript quantification；相比Alignment-based transcript quantification而言，省略了比對這一步驟（有些軟件是只進(jìn)行’輕度’比對），從而直接將reads分配到轉(zhuǎn)錄本上；后者相比前者spliced alignment能極大的減少計算機(jī)資源的使用
現(xiàn)在Salmon支持兩種模式將reads mapping到轉(zhuǎn)錄本上：
第一種是quasi-mapping-based mode，其是一種最新的以及快速準(zhǔn)確的定量轉(zhuǎn)錄本的方法，不需要像傳統(tǒng)的那樣需要通過比對才能將reads mapping到轉(zhuǎn)錄本上；
第二種則是alignment-based mode，跟類似RSEM軟件一樣，提供比對后的bam/sam文件即可對轉(zhuǎn)錄本進(jìn)行定量
第一種方法需要對轉(zhuǎn)錄本建index，第二種方法則不需要
所以也許可以吧！
Salmon到這里下載
下好后放到src文件夾下

tar zxvf salmon-1.1.0_linux_x86_64.tar.gz
cd salmon-latest_linux_x86_64/bin

先試試第一種方法`quasi-mapping-based mode`：

建立索引 2. 對reads進(jìn)行生物學(xué)定量

salmon index -t 你的路徑/hg19.fa -i 你的路徑

index 代表建立索引
-t .fa文件的路徑
-i 索引存放路徑

顯示這樣

然后我突然明白，數(shù)據(jù)需要用轉(zhuǎn)錄本?。。。?br> ……

后來又嘗試了Prokka注釋

但是失敗在了這一步。。。我真的放棄了:-(

2020.6.6更新
在實驗室服務(wù)器上試試

#新建一個index文件存放index
cd ./index/
wget -O hg38.tar.gz https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
tar -zxvf *.tar.gz

一條reads可能會擁有多個valid alignments, 在輸出時，并不會輸出所有的alignments, 而是只輸出-k參數(shù)指定的N個alignments，-k參數(shù)的默認(rèn)值為5。
輸出結(jié)果以SAM格式保存，默認(rèn)輸出到屏幕上，可以通過-S參數(shù)指定輸出文件。

hisat2 -t -x 存放index路徑/index/hg38/genome路徑/genome -1 SRR路徑/SRR3589959_1.fastq -2 SRR路徑/SRR3589959_2.fastq -S 輸出路徑/SRR3589959.sam

-S表示輸出sam文件
只需要寫genome就好，它識別的是前綴，不是某一個文件！
??????可以批量轉(zhuǎn)，像我之前把.fastq都放在一個文件夾里了，并且是雙端測序，所以↓

for i in `seq ？ ？`
do
#將fastq文件轉(zhuǎn)換成sam文件
    hisat2 -t -p 8 -x index的路徑/ctrl/index/hg38/genome -1 sra數(shù)據(jù)的路徑/sra/SRR1234${i}_1.fastq -2 sra數(shù)據(jù)的路徑/sra/SRR1234${i}_2.fastq -S SRR1234${i}.sam
done

這段可以直接在linux命令行運行，或者可以寫成腳本運行
??表示兩個數(shù)字，一個是起始數(shù)字，一個是終止數(shù)字，自己替換

4. Samtools: 轉(zhuǎn)換為`.bam`格式+排序+索引

4.1 `Samtools`將`.sam`文件轉(zhuǎn)為`.bam`文件

for i in `seq ? ?`
do
#將sam文件轉(zhuǎn)換成bam文件
    samtools view -S SRR59618${i}.sam -b > SRR59618${i}.bam
#對bam文件進(jìn)行排序
#剛開始加了參數(shù)-o，運行后電腦終端一直不停跳亂碼的東西，Xshell直接死機(jī)
    samtools sort SRR59618${i}.bam > SRR59618${i}_sorted.bam
#對bam文件建立索引
    samtools index SRR59618${i}_sorted.bam
done

4.2 對`.bam`文件進(jìn)行排序`sort`

samtools sort SRR3589959.bam > SRR3589959_sorted.bam

4.3 對`sorted.bam`文件建立索引

 samtools index SRR3589959_sorted.bam

運行后得到這三種文件

SRR3589959.bam
SRR3589959_sorted.bam
SRR3589959_sorted.bam.bai

#下載hg19_RefSeq.bed文件
$ cd /mnt/f/rna_seq/data/reference/genome/hg19
$ wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz/download
#查看基因組覆蓋率
$ read_distribution.py -i /mnt/f/rna_seq/aligned/SRR3589956.sorted.bam -r /mnt/f/rna_seq/data/reference/genome/hg19/hg19_RefSeq.bed

5. `IGV`可視化結(jié)果

參考這篇

界面簡介

主窗口布局：
1.工具欄tool bar
2.紅色框顯示當(dāng)前顯示的染色體的位置，當(dāng)縮小顯示范圍到整個染色體范圍時，紅色框消失。
3.顯示當(dāng)前查看的染色體序列的長度
4.該窗口顯示測序樣品的測序情況。每一條track代表一個樣品或者一次實驗，顯示的情況包括甲基化、表達(dá)水平、拷貝數(shù)，堿基突變等信息。
5.參考基因組信息
6.track名（即樣品或者實驗名）
7.Attribute names屬性名，即序列信息，如indel、甲基化等。

5.1 選擇`Genomes`→`load genomes from file` 導(dǎo)入下載好的 `.fa`基因組

IGV界面

5.2 載入基因組注釋，選擇`Tools`→`Run igvtools`，進(jìn)入以下`igvtools`窗口：

tools

igvtools

5.3 但是在載入之前需要先將`gff3`進(jìn)行排序，選擇`sort`選項

sort排序

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

RNA-seq摸索：2.sra下載數(shù)據(jù)→fastqc質(zhì)控→hisat2/bowtie2/STAR/salmon比對→Samtools格式轉(zhuǎn)換→IGV可視化結(jié)果

RNA-seq摸索：2.sra下載數(shù)據(jù)→fastqc質(zhì)控→hisat2/bowtie2/STAR/salmon比對→Samtools格式轉(zhuǎn)換→IGV可視化結(jié)果

??前半部分是我的血淚史！正文在中間開始！??

這里都是坑，別學(xué)?。。。≌_的在下一條分界線那里↓

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

1.2 轉(zhuǎn)為`fastq`格式

到這里才完全發(fā)現(xiàn)問題?。。。?！??????

所以就是重頭來過?。。?！??

?所以這里才是正文的開始！！

1 原始數(shù)據(jù)下載

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

注意！?。?！光標(biāo)一定不能在有SRR號的那行，不然會跳過那行??！

1.2 `sra`轉(zhuǎn)成`fastq` (這里要注意不是這么寫?。。?！

因為是雙端測序數(shù)據(jù)?。?！

應(yīng)該加這個參數(shù) `--split-files`

如果需要壓縮，可以加`--gzip`，就會生成`.fastq.gz`文件，解壓縮用`gunzip`

單端測序SE用`-U`參數(shù)

1.3 `md5sum`檢測傳輸?shù)?code>fastq文件是否完整，一般會附帶一個MD5校驗文件

2 質(zhì)控

2.1 `fastqc`生成質(zhì)控報告 fastqc.html和 fastqc.zip格式

2.2 `trim_galore`修剪數(shù)據(jù),用于去除低質(zhì)量和接頭數(shù)據(jù)

3 比對到參考基因組

試試`STAR`

先試試第一種方法`quasi-mapping-based mode`：

4. Samtools: 轉(zhuǎn)換為`.bam`格式+排序+索引

4.1 `Samtools`將`.sam`文件轉(zhuǎn)為`.bam`文件

4.2 對`.bam`文件進(jìn)行排序`sort`

4.3 對`sorted.bam`文件建立索引

運行后得到這三種文件

5. `IGV`可視化結(jié)果

5.1 選擇`Genomes`→`load genomes from file` 導(dǎo)入下載好的 `.fa`基因組

5.2 載入基因組注釋，選擇`Tools`→`Run igvtools`，進(jìn)入以下`igvtools`窗口：

5.3 但是在載入之前需要先將`gff3`進(jìn)行排序，選擇`sort`選項

5.4 輸入注釋文件，點擊`run`，就可以得到`sorted.gff3`文件

5.5 `file`→`load from file`導(dǎo)入`sorted`文件

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

RNA-seq摸索：2.sra下載數(shù)據(jù)→fastqc質(zhì)控→hisat2/bowtie2/STAR/salmon比對→Samtools格式轉(zhuǎn)換→IGV可視化結(jié)果

??前半部分是我的血淚史！正文在中間開始！??

這里都是坑，別學(xué)?。。。≌_的在下一條分界線那里↓

1.1 在SRA數(shù)據(jù)庫下載數(shù)據(jù)

1.2 轉(zhuǎn)為fastq格式

到這里才完全發(fā)現(xiàn)問題?。。。?！??????

所以就是重頭來過?。。?！??

?所以這里才是正文的開始！！

1 原始數(shù)據(jù)下載

1.1 在SRA數(shù)據(jù)庫下載數(shù)據(jù)

注意！?。?！光標(biāo)一定不能在有SRR號的那行，不然會跳過那行??！

1.2 sra轉(zhuǎn)成fastq (這里要注意不是這么寫?。。?！

因為是雙端測序數(shù)據(jù)?。?！

應(yīng)該加這個參數(shù) --split-files

如果需要壓縮，可以加--gzip，就會生成.fastq.gz文件，解壓縮用gunzip

單端測序SE用-U參數(shù)

1.3 md5sum檢測傳輸?shù)?code>fastq文件是否完整，一般會附帶一個MD5校驗文件

2 質(zhì)控

2.1 fastqc生成質(zhì)控報告 fastqc.html和 fastqc.zip格式

2.2 trim_galore修剪數(shù)據(jù),用于去除低質(zhì)量和接頭數(shù)據(jù)

3 比對到參考基因組

試試STAR

先試試第一種方法quasi-mapping-based mode：

4. Samtools: 轉(zhuǎn)換為.bam格式+排序+索引

4.1 Samtools將.sam文件轉(zhuǎn)為.bam文件

4.2 對.bam文件進(jìn)行排序sort

4.3 對sorted.bam文件建立索引

運行后得到這三種文件

5. IGV可視化結(jié)果

5.1 選擇Genomes→load genomes from file 導(dǎo)入下載好的 .fa基因組

5.2 載入基因組注釋，選擇Tools→Run igvtools，進(jìn)入以下igvtools窗口：

5.3 但是在載入之前需要先將gff3進(jìn)行排序，選擇sort選項

5.4 輸入注釋文件，點擊run，就可以得到sorted.gff3文件

5.5 file→load from file導(dǎo)入sorted文件

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

??前半部分是我的血淚史！正文在中間開始！??

這里都是坑，別學(xué)?。。。≌_的在下一條分界線那里↓

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

1.2 轉(zhuǎn)為`fastq`格式

到這里才完全發(fā)現(xiàn)問題?。。。?！??????

所以就是重頭來過?。。?！??

1.1 在`SRA`數(shù)據(jù)庫下載數(shù)據(jù)

注意！?。?！光標(biāo)一定不能在有SRR號的那行，不然會跳過那行??！

1.2 `sra`轉(zhuǎn)成`fastq` (這里要注意不是這么寫?。。?！

因為是雙端測序數(shù)據(jù)?。?！

應(yīng)該加這個參數(shù) `--split-files`

如果需要壓縮，可以加`--gzip`，就會生成`.fastq.gz`文件，解壓縮用`gunzip`

單端測序SE用`-U`參數(shù)

1.3 `md5sum`檢測傳輸?shù)?code>fastq文件是否完整，一般會附帶一個MD5校驗文件

2.1 `fastqc`生成質(zhì)控報告 fastqc.html和 fastqc.zip格式

2.2 `trim_galore`修剪數(shù)據(jù),用于去除低質(zhì)量和接頭數(shù)據(jù)

試試`STAR`

先試試第一種方法`quasi-mapping-based mode`：

4. Samtools: 轉(zhuǎn)換為`.bam`格式+排序+索引

4.1 `Samtools`將`.sam`文件轉(zhuǎn)為`.bam`文件

4.2 對`.bam`文件進(jìn)行排序`sort`

4.3 對`sorted.bam`文件建立索引

5. `IGV`可視化結(jié)果

5.1 選擇`Genomes`→`load genomes from file` 導(dǎo)入下載好的 `.fa`基因組

5.2 載入基因組注釋，選擇`Tools`→`Run igvtools`，進(jìn)入以下`igvtools`窗口：

5.3 但是在載入之前需要先將`gff3`進(jìn)行排序，選擇`sort`選項

5.4 輸入注釋文件，點擊`run`，就可以得到`sorted.gff3`文件

5.5 `file`→`load from file`導(dǎo)入`sorted`文件