1. fastqc 檢查下數(shù)據(jù)質(zhì)量

結(jié)果顯示已去除接頭

2. flash 雙端測(cè)序reads 進(jìn)行拼接

flash -t 6 MTO1_IP.1.fq MTO1_IP.2.fq -p 33 -r 150 -s 100 -o MTO1_IP_merge

生成了6個(gè)文件

image.png

FLASH拼接默認(rèn)輸出6個(gè)結(jié)果文件：

extendeFrags.fastq 為拼接后的擴(kuò)增片段序列文件；

output.flash.log 為日志文件，詳細(xì)記錄了拼接過(guò)程中的參數(shù)和拼接統(tǒng)計(jì)的數(shù)據(jù)；

output.hist 為拼接后的reads長(zhǎng)度的統(tǒng)計(jì)信息文件；

image.png

output.histogram 為拼接后的reads長(zhǎng)度直方圖文件；

image.png

output.notCombined_1.fastq 為拼接不上的reads1序列文件；

output.notCombined_2.fastq 為拼接不上的reads2序列文件；
沒(méi)匹配上的文件根據(jù)需要再進(jìn)行下步分析

3 fastq_quality_filter 去除低質(zhì)量的reads

fastq_quality_filter -q 20 -p 80 -i MTO1_IP_merge.extendedFrags.fastq -o MTO1_IP_mer_clean.fq -Q 33

image.png

4 bowtie2 對(duì)比

bowtie2-build 建立索引, ensemble ftp下載人線粒體序列fasta文件和gff3注釋文件

bowtie2-build ./index/Homo_sapiens.GRCh38.dna.chromosome.MT.fa ./index/human_mt_index

開(kāi)始比對(duì)

bowtie2 -p 10 -x ./index/human_mt_index  MTO1_IP_mer_clean.fq \
| samtools sort -O bam -@ 10 -o - > MTO1_IP_mer_clean.bam

得到bam文件

image.png

5 計(jì)算counts

samtools idxstats適合標(biāo)注的序列進(jìn)行的counts計(jì)數(shù)，如小RNA或者tRNA的序列進(jìn)行索引和匹對(duì)，該結(jié)果得出來(lái)得結(jié)果如下

samtools idxstats MTO1_IP_mer_clean.bam | cut -f 1,3 > MTO1_IP_mer_clean.bam.counts.txt

image.png

對(duì)線粒體16596bp 序列沒(méi)有進(jìn)行注釋。

6 stringtie 計(jì)算count.txt (包括FPKM和RPM)

sudo apt install stringtie

stringtie MTO1_IP_mer_clean.bam -p 16 -G ./gtf/Danio_rerio_Ensemble_97.gtf -B -o ./A/TRANHOM.gtf

這里生成了結(jié)果gtf文件和ballgown需要的.ctab文件，還有基因的表達(dá)量文件gene_abund.tab，該文件包括基因的表達(dá)量FPKM以及TPM等。當(dāng)然如果你想要轉(zhuǎn)錄本的表達(dá)量，直接打開(kāi)t_data.ctab這個(gè)文件，這里面有轉(zhuǎn)錄本的FPKM值。

image.png

image.png

#################################################
多個(gè)轉(zhuǎn)錄本 stringtie --merge 合并多個(gè)gtf文件

image.png

$ stringtie --merge -p 8 -G Danio_rerio_Ensemble_97.gtf -o \
./A/stringtie_merged.gtf ./A/mergelist.txt.txt
或
stringtie --merge -p 8 -G ./GTF/Danio_rerio_Ensemble_97.gtf -o \
./A/stringtie_merged_MCK_HOM.gtf ./a/mergelist.txt

###########################################

7. featureCounts也可以對(duì)bam文件計(jì)算counts數(shù)

featureCounts -T 6 -t exon -g exon_id -a ./gtf/Homo_sapiens.GRCh38.107.chromosome.MT.gff3/Homo_sapiens.GRCh38.107.chromosome.M
T.gff3 -o MTO1_IP.ele.txt MTO1_IP_mer_clean.bam

-g exon_id 需要根據(jù)實(shí)際的gff文件進(jìn)行選擇