转录组测序发现唯一比对率特别低,进一步将未比对到参考基因组的序列提出来比对到NCBI NT库。发现多数序列比对到了Mycoplasma hyorhinis,怀疑是支原体污染,所以我们准备下载Mycoplasma hyorhinis序列参考基因组文件重新比对。
进入NCBI网站,搜索Mycoplasma hyorhinis(https://www.ncbi.nlm.nih.gov/genome/?term=Mycoplasma+hyorhinis),下载gff和fa文件。
下载建索引参考序列文件 gff和fa文件
使用cufflinks里面gffread将gff转为gtf文件。
gffread GCF_000383515.1_ASM38351v1_genomic.gff -T -o \
GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf
使用STAR软件构建索引。STAR建索引特别吃内存,能把你服务器内存用完,然后报错。此时你就要根据你的内存设置limitGenomeGenerateRAM参数,此处设置比所需内存高一点点,如果低于所需内存也会报错。另外线程数可以设置高一点,内存消耗会因为线程数变多而增加,不过不用担心,并不会成倍增加,40个线程内存消耗也就增加了10%。
mkdir S14_INDEX
STAR --runMode genomeGenerate --runThreadN 10 --genomeDir S14_INDEX/ \
--genomeFastaFiles GCF_000383515.1_ASM38351v1_genomic.fna \
--sjdbGTFfile GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf \
--sjdbOverhang 149 #reads长度-1
开始比对
R1=SampleName_R1_001.fastq.gz
outdir=./StepAlign/
gtffile=./star_index/GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf
STAR_index=./star_index/S14_INDEX
R2=${R1//R1_001.paired.fastq.gz/R2_001.paired.fastq.gz}
R1File=`basename ${R1}`
prefix=${R1File//.R1_001.paired.fastq.gz/}
STAR \
--genomeDir ${STAR_index} \
--sjdbGTFfile ${gtffile} \
--limitBAMsortRAM 40000000000 \
--runThreadN 12 \
--limitIObufferSize 500000000 \
--outFilterType BySJout \
--outFilterMismatchNmax 999 \
--outFilterMismatchNoverLmax 0.04 \
--outFilterMultimapNmax 20 \
--outFilterMatchNminOverLread 0.66 \
--outFilterIntronMotifs None \
--outSJfilterReads All \
--outSAMtype BAM SortedByCoordinate \
--outSAMunmapped Within \
--outSAMstrandField intronMotif \
--outSAMattrRGline ID:${prefix} SM:${prefix} PL:Illumina \
--alignSJoverhangMin 8 \
--alignSJDBoverhangMin 1 \
--alignIntronMin 20 \
--alignIntronMax 1000000 \
--alignMatesGapMax 1000000 \
--chimSegmentMin 15 \
--chimJunctionOverhangMin 15 \
--chimScoreMin 0 \
--chimScoreDropMax 20 \
--chimScoreSeparation 10 \
--chimScoreJunctionNonGTAG -1 \
--quantMode TranscriptomeSAM \
--quantTranscriptomeBan IndelSoftclipSingleend \
--outReadsUnmapped Fastx \
--readFilesIn ${R1} ${R2} \
--readFilesCommand zcat \
--outFileNamePrefix ${outdir}/${prefix}.
echo "finish to anign "
网友评论